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SPIS TREŚCI 


Rozdział 1 


Specyfika obliczeń 
numerycznych 


Matematyka obliczeniowa jest dziedziną wiedzy zajmującą się proble- 
mami obliczeniowymi i konstrukcją algorytmów rozwiązywania zadań 
matematycznych. 

Aby w ogóle mówić w problemie obliczeniowym, musimy najpierw 


e określić dane początkowe i cel obliczeń, czyli dokładnie sformuło- 
wać zadanie, 


e określić środki obliczeniowe dzięki którym chcemy osiągnąć cel, 
czyli zdefiniować model obliczeniowy. 


1.1 Zadania numeryczne, przykłady 


Sformułowanie zadania polega na sprecyzowaniu tego, co mamy i co 
chcemy uzyskać. Formalnie polega to na zdefiniowaniu przestrzeni da- 
nych F, przestrzeni wyników G i wskazaniu zależności f > g, gdzie 
feFig€G, między danymi a wynikami. Zależność tę repezentuje 
operator rozwiązania 


S:F—0G. 


Będziemy zainteresowani przede wszystkim zadaniami numerycz- 
nymi. Są to zadania, dla których wynikiem g jest zawsze skończony 
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ciąg liczb rzeczywistych, czyli 
G C URZ 


(przyjmujemy przy tym, że R? = pusty). Zbiór danych może być w 
zasadzie dowolny, ale my będziemy dla uproszczenia rozpatrywać tylko 
zadania dla których F C R”, albo ogólniej, dla których F jest pewną 
klasą funkcji zdefiniowanych na ustalonym zbiorze D C R4. 


Przykład 1.1 (Równanie kwadratowe.) Załóżmy, że chcemy obliczyć 
wszystkie rzeczywiste pierwiastki równania 


m =2ptd= 0, 
dla danych liczb rzeczywistych pi q. Wtedy F = R, G=R"VUR'UR>, 


oraz 
pusty A<0, 
5(p,q) = 4 P Asi (1.1) 
gdzie A = p? — q. 
Przykład 1.2 (Układ równań liniowych.) Rozpatrzmy zadanie roz- 


wiązywania układu równań liniowych 


Aż = b, 
dla nieosobliwej macierzy A = (a, ;);,_, i wektora b= (b;)j_1. Wtedy 
F =4(A,b)eR""": det A40), 


G = R”, oraz z F 
S(A,b) = A'b. 


Przykład 1.3 (Całkowanie.) Dla danej funkcji ciągłej f : [a,b] > R, 
chcemy obliczyć całkę oznaczoną 


W tym przypadku G = R, ale dane stanowi pewna klasa funkcji cią- 
głych określonych na odcinku [a,b|, tzn. F c C(fa,b)). Oczywiście 5 
jest operatorem całkowania, 5 = I. 
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1.2 Model obliczeniowy 


Aby zdefiniować nasz model obliczeniowy, posłużymy się pojęciem pro- 
gramu. Zastosujemy przy tym notację podobną do tej z języka progra- 
mowania Pascal. 


Program składa się z deklaracji, czyli opisu obiektów, których bę- 
dziemy używać, oraz z poleceń (instrukcji), czyli opisu akcji, które bę- 
dziemy wykonywać. 

Dostępnymi obiektami są stałe i zmienne typu całkowitego (inte- 
ger), rzeczywistego (real), albo logicznego (Boolean). Zmienne mogą 
być grupowane w wektory albo tablice. 


Polecenia dzielimy na proste i złożone. Poleceniem prostym jest 


e  PODSTAWIENIE 
ZSZ) 


gdzie z jest zmienną, a W jest wyrażeniem o wartościach tego samego 
typu co z. 


Wyrażeniem jest pojedyncza stała lub zmienna, albo złożenie skoń- 
czonej liczby operacji elementarnych na wyrażeniach. Operacje elemen- 
tarne to: 


arytmetyczno-arytmetyczne: x > —1, (x,y) > x +y, (ay) + 


c—y, (x,y) > w*y, (wy) > z/yy # 0, gdzie x,y są stałymi 
lub zmiennymi liczbowymi, 


arytmetyczno-logiczne: (z,y) > z < y, (x,y) z < y, (x,y) > 
xz = y, (x,y) > x ź y, gdzie x,y są stałymi lub zmiennymi 
liczbowymi, 


logiczno—logiczne: p> notp, (p,q) > pandq, (p,q) > por q, gdzie 
p,q są stałymi lub zmiennymi logicznymi. 


Dla niektórych zadań wygodnie jest (a czasem koniecznie) uzu- 
pełnić ten zbiór o dodatkowe operacje, takie jak obliczanie wartości 
niektórych funkcji standardowych (,/,cos(),sin(), exp(),log(), itp.), 
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czy nawet funkcji bardziej skomplikowanych. Na przykład, zastosowa- 
nie “szkolnych” wzorów na obliczanie pierwiatków równania kwadrato- 
wego z Przykładu 1.1 byłoby niemożliwe, gdyby pierwiastkowanie było 
niemożliwe. Należy jednak przy tym pamiętać, że w praktyce funk- 
cje standardowe (o ile są dopuszczalne) są obliczane używając czterech 
podstawowych operacji arytmetycznych. 

Mamy trzy podstawowe polecenia złożone. 


e WARUNKOWE 
if W then A else Ags; 
gdzie W jest wyrażeniem o wartościach logicznych, a A; i A> są pole- 
ceniami, przy czym dopuszczamy polecenia puste. 


e POWTARZANE 
while W do A; 


gdzie W jest wyrażeniem o wartościach logicznych, a A jest poleceniem. 


e  KOMBINOWANE 
begin A; A>; ... ; An end; 


gdzie A; są poleceniami. 

Na podstawie tych trzech poleceń można tworzyć inne, takie jak 
pętle for i repeat, czy case, itd. 

Mamy też dwa szczególne polecenia, które odpowiadają za ” wejście” 
i wyjście”. 

© WPROWADZANIE DANYCH 


IN (x,t); 


gdzie x jest zmienną rzeczywistą, a t “adresem” pewnego funkcjonału 
L: F — R należącym to pewnego zbioru T. W wyniku wykonania tego 
polecenia w zmiennej z zostaje umieszczona wartość L,(f). 

Polecenie to pozwala zdobyć informację o danej f. Jeśli F = R” to 
zwykle mamy T = {1,2,... n} i L;(f) = fi, co w praktyce odpowiada 
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wczytaniu i-tej współrzędnej wektora danych. W szczególności, ciąg po- 
leceń TN(zfi], i), i = 1,2,...,n, pozwala uzyskać pełną informację o 
f. Jeśli zaś F jest pewną klasą funkcji f : [a,b] — R, to możemy mieć 
np. T = [a,b| i L,(f) = f(t). W tym przypadku, wykonanie polecenia 
IN (x,t) odpowiada w praktyce skorzystaniu ze specjalnej procedury 
(albo urządzenia zewnętrznego) obliczającej (mierzącego) wartość funk- 
cji f w punkcie t. 


e WYPROWADZANIE WYNIKÓW 
OUT (W); 


gdzie W jest wyrażeniem o wartościach rzeczywistych. Polecenie to po- 
zwala “wskazać” kolejną współrzędną wyniku. 


Zakładamy, że na początku procesu obliczeniowego wartości wszyst- 
kich zmiennych są nieokreślone, oraz że dla dowolnych danych wykona- 
nie programu wymaga wykonania skończonej liczby operacji elementar- 
nych. Wynikiem obliczeń jest skońćzony ciąg liczb rzeczywistych (albo 
puste), którego kolejne współrzędne pokazywane są poleceniem OUT. 


1.3 Arytmetyka fi, 


Przedstawiony model obliczeniowy jest modelem idealistycznym, tzn. 
zakłada on, że wszystkie operacje są wykonywane bezbłędnie. Dlatego 
w tym przypadku będziemy mówić o arytmetyce idealnej. W praktyce 
jednak, np. wykonując obliczenia na maszynie cyfrowej, operacje aryt- 
metyczne na liczbach rzeczywistych wykonywane są z pewnym błędem. 
Matematycznym modelem arytmetyki maszyny cyfrowej jest arytme- 
tyka fi, (albo arytmetyka zmiennoprzecinkowa), którą teraz przedsta- 
wimy. 


Dowolną liczbę rzeczywistą z Æ 0 można przedstawić w postaci 
x = s- 10°- m, 


gdzie s € {—1, 1} jest znakiem, c € Z cechą, am € [0.1,1.0) mantysą 
liczby xz. Zauważmy, że taki rozkład jest jednoznaczny i odpowiada prze- 
suwaniu przecinka w rozwinięciu dziesiętnym liczby do pierwszej cyfry 
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znaczącej, tj. różnej od zera. Mantysa ma w ogólności nieskończenie 
wiele cyfr c; w swoim rozwinięciu dziesiętnym, 


m = X g1, 
j=1 
gdzie c; € {0,1,2,...,9} ic, 4 0. Jako taka nie może więc być zapa- 
miętana dokładnie w pamięci maszyny cyfrowej. Zakładając, że pamię- 
tamy tylko t cyfr znaczących i ostatnią cyfrę zaokrąglamy, zamiast m 
pamiętamy 
t 
m = X G10 + dzl0””, 
j=1 

gdzie 

= 0 gdy 0 < cı < 4}, 

ace" eh gdy 5 < ci < 9}. 


Stąd z będzie reprezentowana przez liczbę 
rd(z) = s- 10°- m. 
Między liczbą z (x # 0) a jej reprezentacją rd(x) zachodzi więc 
nierówność 
|x — rd(z)|  |]m— m] za: 10-69 
|z| Spa 107 


= 5.107, 


charakteryzująca dokładność arytmetyki. Liczby rzeczywiste są więc w 
arytmetyce fl, reprezentowane z dokładnością względną v = 5:107*. (We 
współczesnych maszynach cyfrowych w pojedynczej precyzji mamy co 
najmniej t = 8.) Zauważmy, że t jest największą spośród liczb natural- 
nych j spełniających bierówność rd(1 + 107) > 1. 

Ostatnią nierówność wygodnie jest zapisać w równoważny sposób 
jako 

rd(x) = z(1 +€), gdzie |e| < v. 


W arytmetyce fl, zakładamy ponadto, że działania arytmetyczne na 
liczbach rzeczywistych (a raczej na ich reprezentacjach) są wykonywane 
dokładnie i tylko wynik jest zaokrąglany. Mamy więc 


A,(zDy) = rd ( rd(x) O rd(y) ), 


1.3. ARYTMETYKA FL, T 


gdzie O € (+, —,*, /}, Ogólniej, jeśli Wı i W2 są wyrażeniami o war- 
tościach rzeczywistych, to dla dowolnych wartości zmiennych 


f,(WEW2) = rd (A, (W1) DA, W2)). 


Zwykle dla prostoty będziemy również zakładać podobną zależność dla 
niektórych funkcji standardowych, o ile należą one do zbioru operacji 
elementarnych (chociaż w rzeczywistości są one obliczane przez pro- 
cedury używające czterech podstawowych operacji arytmetycznych). I 
tak będziemy mieć np. 


A(W) = (VAOW)A+ a). 
AW) = (eos(ft,(W)) (+ 2), 


gdzie |e;| < v, oraz 8, < K,v i K; są “niewielkimi” stąłymi. 

Podobnie, jeśli A jest operatorem porównania, A € {<, <, =, Æ}, 
to wartością wyrażenia logicznego W1 AW., w fi, jest dokładna wartość 
wyrażenia fi,(W)Af,(W>). 


Uwagi i uzupełnienia 


U. 1.1 W maszynie cyfrowej cecha c liczby rzeczywistej nie może oczywi- 
ście mieć dowolnie dużej wartości bezwzględnej, |c| < Cmax. Powoduje to 
powstanie zjawiska nadmiaru gdy c > Cmax, oraz zjawiska niedomiaru gdy 
C < —Cpmax. W pierwszym przypadku liczba jest tak duża (co do modułu), 
że nie zawiera się w przedziale liczb reprezentowalnych, a w drugim jest tak 
mała, że musi być reprentowana przez zero, przy czym błąd względny repre- 
zentacji wynosi wtedy 1 a nie v. W dalszych rozważaniach zjawiska nadmiaru 
i niedomiaru będziemy dla uproszczenia zaniedbywać. 


U. 1.2 W następnych rozdziałach często będziemy się posługiwać normami 
wektorów £ = (u;);_, € R” i macierzy A = (a,.;);,—, € R"*", a w szcze- 
gólności błędami reprezentacji wektorów i macierzy w normie. Najczęściej 
używanymi normami wektorowymi będą normy p-te, 


A 1/p 
lzi = Iž = (Żar) ,  1śp<+o, 
j=1 
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Oraz 


Ilo = „im [Flp = max lezl 


Normą macierzową jest norma euklidesowa (albo Frobeniusa) 


a także normy indukowane przez normy wektorowe (np. przez normy p-te) 


Az] s 
—— = sup |Az||. 
IE lel=1 


IA]| = sup 
x#0 


Jeśli norma macierzowa jest indukowana przez normę wektorową, to dla 
dowolnego wektora mamy 


IAZ < AIZI] 


U. 1.3 Przypomnijmy, że w przestrzeniach liniowych skończenie wymiaro- 
wych (a więc także w R” i w przestrzeni macierzy wymiaru n x n) każde 
dwie normy są równoważne. To znaczy, że jeśli mamy dwie normy ||- || i J|- 
w przestrzeni skończenie wymiarowej X, to istnieją stałe 0 < Kı < Kə < œo 
takie, że 

Kı |z| < æl < Kaja,  Vz€X. 


W szczególności dla x € R” mamy 


lžlo < lin < nll, (1.2) 
lžllo < |F]2 < va lllo, (1.3) 
1 
ym Ml < lže < ||7||1, (1.4) 
a dla A = (ai j) j=1 mamy 
lAl < IA] l2 < Als < v” lAl, (1.5) 


gdzie |A| = (laig i= 
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Cwiczenia 
Ćw. 1.1 Podać przykłady liczb z i y, które są dokładnie reprezentowane w 
fir, ale filz y) TUY, gdzie E {+, =, *, /} 


Ćw. 1.2 Uzasadnić, że jeśli £ < y to rd(x) < rd(y) oraz fi,(z/y) < 1. 


Podać przykład wyrażenia W takiego, że dla pewnych wartości zmiennych 
W <0, ale f (W) > 0. 


Ćw. 1.3 Jak zabezpieczyć się przed powstaniem nadmiaru albo niedomiaru 
przy obliczaniu wyrażenia Va? + y?, gdy x i y leżą w przedziale liczb repre- 
zentowalnych w fl,, ale z? lub y? nie? 


Ćw. 1.4 Uzasadnić nierówności (1.2)-(1.5). 


Ćw. 1.5 Pokazać, że dla macierzy A = (ai), 1 € R” mamy 


n 
Allo = ESPON 


Oraz 


n 
All = IA" = max X |a|. 
i<j<n =? 


Ćw. 1.6 Dla wektora 7 = (c;)j_q E R”, niech rd(¥) = (rd(z;));_4. Poka- 
zać, że 
Z- rd(2)||p < v|£]|p 


dla 1 < p < œ. 


Ćw. 1.7 Dla macierzy A = (aij)ęj=1, niech rd(A) = (rd(a;.;))ę,=1. Poka- 
zać, że 


[A= ra(A)|p < v ||Allp, 


dla p = 1, œ, oraz 


IA 


[.4— rd(A)||2 < |A- ra(A)lle < vl|Alle < vnr ||All2. 
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Rozdział 2 
Algorytmy i ich własności 


Algorytm to dokładnie określona i dozwolona w danym modelu oblicze- 
niowym sekwencja akcji, pozwalająca na rozwiązanie danego zadania 
(w sposób dokładny lub przybliżony). 


2.1 Rozkład algorytmu względem infor- 
macji 
Z każdym algorytmem związany jest operator 
$: F— G, 


taki że $(f) jest wynikiem działania algorytmu w arytmetyce idealnej 
dla danej f. 

Zauważmy, że wynik 9(f) działania algorytmu nie zależy bezpo- 
średnio od f, ale raczej od informacji o f (uzyskanej dzięki poleceniu 
TN). Informacja ta może być pełna albo tylko częściowa. Informacja 
jest pełna gdy, np. f=(f1,..., fn) € R” i wczytamy wszystkie współ- 
rzędne f,. Informacja może być częściowa, gdy f jest funkcją. Wtedy 
wiele danych może posiadać tę samą informację, co łatwo zaobserwować 
na przykładzie zadania całkowania. 

Niech N : F — UX oR” będzie operatorem informacji, tzn. 


N(F) = (Yy1,Y2,---;Yn) 


11 
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jest informacją o f zebraną przy idealnej realizacji algorytmu. Za- 
uważmy, że nformacja jest pełna gdy N jest przekształceniem różnowar- 
tościowym, tzn. jeśli fi # f2 implikuje N(f1) Æ N(f2). W przeciwnym 
przypadku mamy do czynienia z informacją częściową. 

Każdy algorytm © może być przedstawiony jako złożenie operatora 
informacji i pewnego operatora y : N(F) — G zdefiniowanego równo- 
ścią 

p(N(F)) = *(7). 

Zauważmy, że w przypadku informacji częściowej zwykle nie istnieje 
algorytm dający dokładne rozwiązanie zadania dla każdej danej f € F, 
ponieważ dla danych o tej samej informacji mogą istnieć różne rozwią- 
zania. 


2.2 Problem wyboru algorytmu 


Wybór algorytmu jest najistotniejszą częścią całego procesu numerycz- 
nego rozwiązywania zadania. Kierujemy się przy tym przede wszystkim 
następującymi kryteriami: 


e dokładnością algorytmu, 
e złożonością algorytmu, 
e własnościami numerycznymi algorytmu. 


Przez dokładność algorytmu rozumiemy różnicę między rozwiąza- 
niem dokładnym S(f), a rozwiązaniem ©(f) dawanym przez algorytm 
w arytmetyce idealnej. Jeśli *(f) = S(f), Vf € F, to algorytm nazy- 
wamy dokładnym. 

Mówiąc o złożoności, mamy na myśli złożoność pamięciową (zwy- 
kle jest to liczba stałych i zmiennych używanych przez algorytm), jak 
również złożoność obliczeniową. Na złożoność obliczeniową algorytmu 
dla danej f składa się koszt uzyskania infomacji y = N(f) (zwykle jest 
on proporcjonalny do liczby wywołań polecenia ZN ), oraz koszt kombi- 
natoryczny przetworzenia tej informacji, aż do uzyskania wyniku p(y). 
Koszt kombinatoryczny zwykle mierzymy liczbą operacji arytmetycz- 
nych wykonywanych przez algorytm. 
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Przez własności numeryczne algorytmu rozumiemy jego własności 
przy realizacji w arytmetyce fl„. Temu ważnemu tematowi poświęcimy 
teraz osobny paragraf. 


2.3 Numeryczna poprawność algorytmu 


Pożądane jest, aby algorytm dawał “dobry” wynik zarówno w arytme- 
tyce idealnej, jak i w arytmetyce f/,. Niestety, jak zobaczymy, nie zawsze 
jest to możliwe. Nawet jeśli algorytm jest dokładny to w wyniku jego 
realizacji w fl, możemy otrzymać wynik fl,(®(f)) daleko odbiegający 
od S(f). W szczególności, prawie zawsze mamy 


SUV RAEI 


Zauważmy również, że o ile z reguły znamy dokładne zachowanie się al- 
gorytmu w arytmetyce idealnej dla danej informacji, to nie można tego 
samego powiedzieć o jego zachowaniu się w arytmetyce fl,. W związku 
z tym powstaje pytanie, jak kontrolować błąd algorytmów wynikający z 
błędów zaokrągleń i jakie algorytmy uznamy za te o najwyższej jakości 
numerycznej. 


Istnienie błędów reprezentacji liczb rzeczywistych powoduje, że in- 
formacja y = N(f) o danej f nie jest w ogólności reprezentowana do- 
kładnie. Znaczy to, że zamiast na informacji dokładnej, dowolny al- 
gorytm będzie operować na informacji nieco zaburzonej yy, tzn. zabu- 
rzonej na poziomie błędu reprezentacji. Tak samo wynik dawany przez 
algorytm będzie w ogólności zaburzony na poziomie błędu reprezenta- 
cji. W najlepszym więc wypadku wynikiem działania algorytmu w fl, 
będzie (o(y„))„ zamiast p(y). Algorytmy dające tego rodzaju wyniki 
uznamy za posiadające najlepsze własności numeryczne w arytmetyce 
fl, i nazwiemy numerycznie poprawnymi. 

Dokładniej, powiemy, że ciąg rzeczywisty a, = (Gy4,.-.,dyn) (A 
właściwie rodzina ciągów fa„|„) jest nieco zaburzonym ciągiem a = 
(a1,...,aqn), jeśli istnieje stała K taka, że dla wszystkich dostatecznie 
małych v zachodzi 


avj; — a| < K vlajl, 1<jJ<n, 2.1 
j j j 


14 ROZDZIAŁ 2. ALGORYTMY I ICH WŁASNOŚCI 


albo ogólniej 
las — a|| < Kv ljal, (2.2) 
gdzie ||- || jest pewną normą w R”. W pierwszym przypadku mówimy 
o zaburzeniu “po współrzędnych”, a w drugim o zaburzeniu w normie 
I-II. 
Zauważmy, że niewielkie zaburzenia po współrzędnych pociągają za 
sobą niewielkie zaburzenia w normie. Rzeczywiście, jeśli (2.1) to również 


lav — allo = max |avj — aj] < Kv max |a;| = Kv||a||e, 


i korzystając z faktu, że w przestrzeni skończenie wymiarowej wszystkie 
normy są równoważne otrzymujemy dla pewnych stałych K; i Ka 


a, — a]| < Killa — allo < KiK v llall < KKK v ljal, 


czyli nierówność (2.2) ze stałą KoK;K zamiast K. 


Definicja 2.1 Algorytm © rozwiązywania zadania nazywamy nume- 
rycznie poprawnym w zbiorze danych Fo C F, jeśli dla każdej danej 
f € Fo wynik fi,(P(f)) działania algorytmu w arytmetyce fi, można 
zinterpretować jako nieco zaburzony wynik algorytmu w arytmetyce 
idealnej dla nieco zaburzonej informacji y, = (N(f)), E€ N(F) o f, 
przy czym poziom zaburzeń nie zależy od f. 

Formalnie znaczy to, że istnieją stałe K1, Ka, oraz w > 0 takie, że 
spełniony jest następujący warunek. Dla dowolnej v < vo oraz informa- 
cji ye N(Fy) można dobrać y, € N(F') oraz (p(y„)), takie, że 


ly = yll < Kiellyll, 
Ilun), — e)l £ Kar ||ply)||, 


oraz 


fs (D) = As ONF) = e). 


Zauważmy,że jeśli f € R”, N(f) = (fi,.--, fn), oraz algorytm jest 
dokładny, ® = y = S$, to numeryczną poprawność algorytmu można 
równoważnie zapisać jako 


fs (E = (S), 
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2.4 Rola uwarunkowania zadania 


Niech $(.) = o(N(-)) będzie algorytmem numerycznie poprawnym dla 
danych Fo C F. Wtedy jego błąd w fi, można oszacować następująco: 


ISC) — A DIL = IFK) = (ew), I 
WI le) — e) + let) — (e(y)), I 
DIL + 260) — ew) + Kar lleu) 

ISM) -=I + A + Kov) — ew) + Kar liel), 
przy czym |y, — yl| < Kır |lyll. Stąd w szczególności wynika, że jeśli 
algorytm jest numerycznie poprawny i ciągły ze względu na informację 
y, to 


lim [S = ALDI = IS) — Sl 


To znaczy, że dla dostatecznie silnej arytmetyki algorytm będzie się 
zachowywał w fl, prawie tak jak w arytmetyce idealnej. 


Z powyższych wzorów wynika, że błąd w fl, algorytmu numerycznie 
poprawnego zależy w dużym stopniu od: 


e dokładności algorytmu w arytmetyce idealnej, 
e dokładności v arytmetyki fł,, 
e wrażliwości algorytmu na małe względne zaburzenia informacji y. 


Ponieważ dwa pierwsze punkty są raczej oczywiste, poświęcimy tro- 
chę więcej uwagi jedynie trzeciemu. 
Jeśli p spełnia warunek Lipschitza ze stałą L, a dokładniej 


ley) — pl < Lily — yll, 


to 


ISE) — ALEI 
< ISA) -ENI + A+ Kr) Lljy, — yll + Karlil) 
< ISA -EAI + A+ Kv) LKywjjy|| + Kzrllp(y)||. 


W tym przypadku błędy zaokrągleń zwiększają błąd bezwzględny al- 
gorytmu proporcjonalnie do v. 
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Bardziej jednak interesuje nas błąd względny. Wybierzmy “małe” 
n > 0 i przypuśćmy, że 


lo) — p(yji < MKiv max(n, |p(y)|l); 


dla pewnej M niezależnej od y, tzn. błąd względny informscji, ||y, — 
yl| < Kąwllyl|, przenosi się na błąd względny wyniku (w arytmetyce 
idealnej) ze “współczynnikiem wzmocnienia” M, albo na błąd bez- 
względny ze współczynnikiem My. (Zauważmy, że gdybyśmy wzięli 
n = 0 to dla y takiej, że p(y) = 0 musiałoby być p(y,) = 0, co zwykle, 
choć nie zawsze, nie jest prawdą.) Wtedy 


ISC) — ALEI 


IS) — S(ÓI+ + Kay) M Kv max(n, LGD + Korlo) 
IS) - (0 + v (MKA + Kav) + K2) max(n, lel). 


IA 


W szczególności, gdy algorytm jest dokładny i korzysta z pełnej infor- 
macji o f, tzn. S = $ = ọ, to błąd 


ISA — ALCI > 
max, ISC < (MK;(1+ Kav) + Kə) v ~ (MK, + Ky)v. 
Stąd wynika, że jeśli (MK, + Ką)v < 1 to błąd względny algorytmu w 
fi, jest mały, o ile ||S(f)|| > n. Błąd względny jest proporcjonalny do 
dokładności v, arytmetyki fl,, współczynników proporcjonalności K; 
algorytmu numerycznie poprawnego, oraz do wrażliwości M zadania 5 
na małe względne zaburzenia danych. 


Zwróćmy uwagę na istotny fakt, że interesują nas właściwie tylko 
te zaburzenia danych (informacji), które powstają przy analizie algo- 
rytmu numerycznie poprawnego. I tak, jeśli algorytm jest numerycznie 
poprawny z pozornymi zaburzeniami danych w normie, to trzeba zba- 
dać wrażliwość zadania ze względu na zaburzenia danych w normie. 
Jeśli zaś mamy pozorne zaburzenia *po współrzędnych” (co oczywiście 
implikuje pozorne zaburzenia w normie) to wystarczy zbadać uwarun- 
kowanie zadania ze względu na zaburzenia “po współrzędnych”, itd. 

Zadania, które nie są zbyt wrażliwe na “małe” względne zaburze- 
nia danych, tzn. dla których M jest “niewielkie”, nazywamy ogólnie 
zadaniami dobrze uwarunkowanymi. 
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2.5 Przykłady 


Podamy teraz proste przykłady zadań, które mogą być rozwiązane al- 
gorytmami mumerycznie poprawnymi. Oszacujemy też błędy tych al- 
gorytmów. 


2.5.1 Iloczyn skalarny 


Załóżmy. że dla danych ciągów rzeczywistych o ustalonej długości n, 
aj, bj, 1 £ j < n, chcemy obliczyć 


Rozpatrzymy algorytm dokładny zdefiniowany powyższym wzorem i 
korzystający z pełnej informacji o kolejnych współrzednych. 

Oznaczmy przez d, i b, reprezentacje liczb a, i b; w fi,, A; = a;(1+ 
aj), b; = b;(1 + 85), oraz przez 4, i 0, błędy względne powstałe przy 
kolejnych mnożeniach i dodawaniach. Oczywiście |a;|, [8;|, lysh [05] < v. 
Otrzymujemy 


n (aż = (MÓGŁ) + GO +m) J(+) = o 


= (o (abia +q) +b +2) +6) 


Betin} 11) (1+6,) 
= ãba(1 +7)(1 + 82): (1+ ôn) 


= $ ab;(1+e;), 
j=1 


gdzie w przybliżeniu (tzn. gdy v — 0) mamy ļe1| < (n +2) i ļe;| £ (n— 
j+4)v,2 < j < n. Algorytm naturalny jest więc numerycznie poprawny 
w całym zbiorze danych, gdyż wynik otrzymany w fl, można zinterpre- 
tować jako dokładny wynik dla danych a,, = a, ib,, = b;(1 +e;), przy 
czym ||b, — bļlp < (n + 2)v||b||p. 
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Zobaczmy teraz, jak błąd we współrzędnych b; wpływa na błąd 
wyniku. Mamy 


Dab; - AO a;b,)| = | ab - Kabi (1+ ej) | 
= |D eo, 5 5 le;llazi 
j=1 j=1 


(n + 2)v X |a;b;l. 
j=l 


IAR 


Stąd dlay> 0 


| 2% —1 jbj "JM a;b,)| 
AO | X; ajb;|) 


< K,(n+2)v, 


gdzie 
21 lazb;l 


K, = K b) = ; 
E E E) 


Zauważmy, że jeśli iloczyny ajb; są wszystkie dodatnie albo wszyst- 
kie ujemne, to K, = 1, tzn. zadanie jest dobrze uwarunkowane, a błąd 
względny jest zawsze na poziomie co najwyżej nv. W tym przypadku 
algorytm zachowuje się bardzo dobrze, o ile liczba n składników nie 
jest horendalnie duża. W ogólności jednak K,„ może być liczbą dowol- 
nie dużą i wtedy nie możemy być pewni uzyskania dobrego wyniku w 


A, 


2.5.2 Całkowanie 


Zadanie całkowania z Przykładu 1.3 często rozwiązuje się (a raczej przy- 
bliża) stosując formułę 


BH) = ols) = E c 


gdzie informacja y; = f(t;). tj są tu ustalonymi punktami z [a,b], a c; 
ustalonymi współczynnikami rzeczywistymi, niezależnymi od f. Algo- 
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rytmy korzystające z takich formuł nazywamy kwadraturami, a przy- 
kładem jest zwykła suma Riemanna, 


Pę(f) = L 


O; 


1 


I 


gdzie t; należy do przedziału [a + (j — 1)(b — a)/n,a + j(b — a)/n], 
I<j<n. 

Z analizy algorytmu obliczającego iloczyn skalarny wynika, że kwa- 
dratura jest algorytmem numerycznie poprawnym. Mamy bowiem 


fy (2 aiw] = >, czy; + ej), 
j=l j=1 
gdzie |e;| < (n + 1)v. (W porównaniu z iloczynem skalarnym mamy tu 
n+ 1 zamiast n + 2, bo c; nie jest daną i nie podlega zaburzeniu.) Stąd 
błąd bezwzględny kwadratury w fl, można oszacować następująco: 


Uwagi i uzupełnienia 


U. 2.1 Wprowadzimy teraz formalną definicję zastosowanej już wcześniej 
przybliżonej nierówności <. Dla dwóch funkcji piszemy 


[h4(7)| < |h2(7)| 


gdy 
lim sup [ka (7)|/|ka(v)| <1. 


Na przykład, jeśli |h(v)| < X= Kul (Kı > 0), to |h(v)| £ Kiv, co już 
zauważyliśmy w przykładzie z iloczynem skalarnym. Zapis ten wyraża prosty 
fakt, że dla praktycznych wartości v (v < 1078) wyrazy typu v?, v3, itd. są 
tak małe w porównaniu z v, że można je zaniedbać. 
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Łatwo sprawdzić, że przy notacji < zachodzą następujące fakty. Jeśli 
je] < Kv, |eq| < Kiv i |eg| < Kov (gdzie €, £1, €2 są funkcjami v), to 


(1+e1)(1+€2) = (1+ ô), gdzie — |ôļ|< (Kı + Kə)v, (2.3) 


(1+e)1=(1+6), gdzie  |ô| £ Kv, (2.4) 
si 
(+e) = (1+8), gdzie  |ó|< „Ku, (2.5) 
i ogólnie 
(+e) = (1+6), gdzie || < |plv. (2.6) 


U. 2.2 Rozpatrzymy teraz zadanie obliczenia wszystkich pierwiastków rze- 
czywistych równania kwadratowego z Przykładu 1.1. Będziemy zakładać, że 
model obliczeniowy dopuszcza obliczanie pierwiastków kwadratowych z liczb 
nieujemnych oraz fl (yx) = rd(yrd(z)). 

Okazuje się, że nie umiemy pokazać numerycznej poprawności “szkol- 
nego” algorytmu obliczającego pierwiastki równania bezpośrednio ze wzorów 
(1.1). Można jednak pokazać numeryczną poprawność drobnej jego modyfi- 
kacji wykorzystującej wzory Viete'a. 


A :=p*p—gq; 
if (A = 0) then OUT(p) else 
if (A>0) then 


begin 
A1 := sqrt(d); 
if (p > 0) then 
begin 
xl := p + Al; 
z2 := q/zxl; 
end else 
begin 
r2 := p — Al; 
xl := q/x2; 
end; 
OUT (x1); OUT (x2) 
end. 


Mamy bowiem 


f(Alp a) = (PA+0)7(1+e) -a +8))(1+ 22) 
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(1+ 8) 
(r TQ + a)2(1 + e1) 


(P- q(1+8))1+9) = Ap,q(1+5))(1+7), 


) ERN E EER) 


gdzie |ó|, |y| < 4v. Wyróżnik obliczony w fl, jest więc nieco zaburzonym wy- 
różnikiem dokładnym dla danych p i qy = q(1 + ô). W szczególności 


sgn(fi,(A(p,q))) = sgn(A(p,q,)). 


Jeśli p > 0 to 


Gi) = (PO +a) + YA(AG,0)A +e3))(1 + ex) 


(P +a) + y Alp, a) +70 + e3))(1 + e4) 
= r yA, a CEE | (4601-40 


= A(p, 1) (L+), 


gdzie |ej| £ 6v. Zauważmy, że ostatnia równość zachodzi dlatego, że doda- 
jemy liczby tego samego znaku. (Inaczej |e;| mogłaby być dowolnie duża i 
tak byłoby w algorytmie szkolnym.) Dla drugiego pierwiastka mamy 


q(1 + 8) 
f„(c1(p,q)) 


du 


A,(22(p,q)) = f,(z1(p,q)) 


(1 +e5) = (1 + e2), 


gdzie |e2| < 8v. 

Podobny wynik otrzymalibyśmy dla p < 0. Algorytm zmodyfikowany 
jest więc numerycznie poprawny, gdyż otrzymane w fl„ pierwiastki są nieco 
zaburzonymi dokładnymi pierwiatkami dla danych py = p i qv = q(1 + ô). 


Aby oszacować błąd algorytmu, wystarczy zbadać uwarunkowanie zada- 
nia ze względu na zaburzenie danej q, ponieważ pokazaliśmy, że zaburzenia 
p można przenieść na zaburzenia q i wyniku. Niestety, choć algorytm jest nu- 
merycznie poprawny, zaburzenia q mogą sprawić, że nawet znak wyróżnika 
A może być obliczony nieprawidłowo. Na przykład dla p = 1 i q = 1 + 10+: 
mamy A(p,q) = F10t+1, ale A(rd(p), rd(q)) = A(1,1) = 0. Ogólnie 


IA (A(p,4)) — A(p,q)| £ 4u(p” + 2/q|). 
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a więc tylko dla |A(p,q)| = |p? — q| > 4v(p? + 2ļq|) możemy być pewni 
obliczenia właściwego znaku A. Przy tym warunku oraz A > 0 błąd danych 
przenosi się w normie euklidesowej na błąd wyniku następująco: 


(e1(p.9 — 21(p,0)? + (a2(p.9) — «2(p,0)7) 


_ v2|óg| 2 zy2y II 


VP —q+vyP" —qG V? -q 


3 lal / p? 
a v1- q/p max(n/|p|, V2(1 + (1 — q/p*))) 


- max(n, (£1 (p, q)? + x2(p, 0>) 1P). 


Stąd widać, że zadanie jest dobrze uwarunkowane dla q/p? « 1 i może być 
źle uwarunkowane dla q/p? z 1. W ostatnim przypadku nie możemy być 
pewni otrzymania dobrego wyniku w ff,. 


Cwiczenia 
Ćw. 2.1 Pokazać równości (2.3)-(2.6). 
Ćw. 2.2 Niech 0 < aj < a2 < ---aņn. Czy z punktu widzenia błędów w fa 


lepiej jest policzyć sumę tych liczb w kolejności od najmniejszej do najwięk- 
szej czy odwrotnie? 


Ćw. 2.3 Aby obliczyć S (a,b) = a? — b? można zastosować dwa algorytmy: 
P;(a,b) = axa—bx*b oraz $o(a,b) = (a +b) * (a — b). Pokazać, że oba 
algorytmy są numerycznie poprawne, ale drugi z nich wywołuje mniejszy 
błąd względny wyniku w przypadku, gdy rd(a) = ai rd(b) = b. 


Ćw. 2.4 Pokazać, że naturalny algorytm obliczania cosinusa kąta między 
dwoma wektorami d,b € R”, 


Q gd ajbj 


CRIER 


jest numerycznie poprawny. Oszacować błąd względny wyniku w ff,. 


cos(a,b) = 
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Ćw. 2.5 Pokazać, że naturalny algorytm obliczania ||Aż||2 dla danej ma- 
cierzy A € R”? i wektora i € R” jest numerycznie poprawny. Dokładniej, 


AZ(A2]|2) = (A+ EJF, 
gdzie |F||> < 2(n + 2)ynv||Al|2. Ponadto, jeśli A jest nieosobliwa to 


IA,(ILAZ[|2) — AZ| < 2(n + 2vnv ([IAJ2lLA7*|l2) Až 


Ćw. 2.6 Niech © będzie algorytmem numerycznie poprawnym w zbiorze 
danych f € Fo, przy czym dla małych v, fl (®(f)) = o(y), gdzie |y,—y|| < 
Kv||jy|| i K nie zależy od vi f (y= N(f)). Pokazać, że w ogólności © nie 
musi być “numerycznie poprawny po współrzędnych”, tzn. w ogólności nie 
istnieje bezwzględna stała Kı taka, że dla małych v i dla dowolnej f € Fo 


lwi- yl < kvlyjj,  1<j<n, 
gdzie y = (y1, ---, Yn). 


æ E 
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ROZDZIAŁ 2. ALGORYTMY I ICH WŁASNOŚCI 


Rozdział 3 


Układy równań liniowych 


Rozpoczynamy analizę metod dokładnych rozwiązywania układów rów- 
nań liniowych postaci 
AD = b, (3.1) 


gdzie A = (a;;) jest macierzą nieosobliwą (det A # 0) wymiaru 
n x n, ab = (b); jest wektorem z R”. Zakładamy, że informacja 
o zadaniu dana jest przez współczynniki a;; macierzy i współrzędne 
b, wektora. W wyniku powinniśmy uzyskać współrzędne 1; wektora 


. . dk x \N 
rozwiązania £ = (z7), 


n 
ij=1 


A N 


Przypomnijmy, że nieosobliwość macierzy A zapewnia, że rozwiązanie 
z* istnieje i jest wyznaczone jednoznacznie. 


3.1 Układy z macierzą trójkątną 


Szczególnym przypadkiem (3.1) jest układ z macierzą trójkątną A. 
Będą nas szczególnie interesować macierze trójkątne górne, dla których 
aij = 0 gdy i > j, oraz macierze trójkątne dolne z jedynkami na prze- 
kątnej, tzn. a,, = 0, i < j, oraz a;i; = 1. Macierze pierwszego rodzaju 
będziemy oznaczać przez R, a drugiego rodzaju przez L. 

Układ z macierzą trójkątną górną 


Ri=e (3.2) 
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R= (rij), C= (cj), można rozwiązać stosując algorytm: 


Th = Cn/Tna; 

for i := n—1 downto 1 do 
r n 

T} = (ci = aa rajat) [Tii 


(Algorytm ten jest wykonalny, ponieważ nieosobliwość macierzy impli- 
kuje, że r;, Æ 0, Vi.) Podobnie, układ Li = č rozwiązujemy algoryt- 
mem: 


ti =Lch 

for i :=2 to n do 
a a 1—1 * 

Ti = G — X= lijt}. 


Oba algorytmy wymagają rzędu n?/2 mnożeń lub dzieleń i n?/2 doda- 
wań lub odejmowań, a więc rzędu n? wszystkich działań arytmetycz- 
nych. 


3.2 Eliminacja Gaussa 


Eliminacja Gaussa jest algorytmem dokładnym rozwiązywania ukła- 
dów równań z macierzą pełną Á, polegającym na sprowadzeniu układu 
wyjściowego (3.1) do układu trójkątnego (3.2). 

W postaci rozwiniętej układ (3.1) ma postać 


d1iTi T Grodo T *': T dlntn = bı 
d21%1 T dod T ° T dąntln = ba 
dn 1X1 T dn2l2 T ° T Annin = bn. 


Zauważmy, że jeśli od i-tego wiersza tego układu odejmiemy wiersz 
pierwszy pomnożony przez l1 = Q1/011, dla i = 2,38,...,n (zakła- 
damy, że a, Æ 0) to otrzymamy układ równoważny A® z = bl), 


ayiti + Qoto + © + Ainín = DI 


a e e a a 5 M) 


aaa +o + aea = LP, 
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z wyrugowaną zmienną xı z równań od drugiego do n-tego. Teraz z 
kolei możemy od i-tego wiersza odjąć wiersz drugi pomnożony przez 
li2 = af; land, dla i = 3,4,...,n (o ile as? Æ 0), aby otrzymać układ 
ABZ = bO) z wyrugowanymi zmiennymi 2, i £2 z równań od trzeciego 
do n-tego. Postępując tak dalej otrzymujemy po n — 1 krokach układ 
AG-Dz = b") postaci 


OWCA + aj 412 o ew k 26 = PI 
ała + + cdi, = b 
atr Den = O=), 


Układ ten jest, trojkatny górny, RZ = €, z macierzą R = A") i 
wektorem £ = bl” , można go więc rozwiązać znaną już metodą. 


Kolejne operacje w eliminacji Gaussa można zapisać następująco: 


for k:=1 to n—1l do 
for 1: k+l to n do 


begin 
ię = aję" Jagy 
DB = „ED aka lik x W 
for j:=k+1 to n do 
k k—1 k—1 
Gy = Są 2 li EFG), ) 
end. 


Zauważmy, że k-ty krok wymaga rzędu (n — k)? dodawań lub odejmo- 
wań i BĘ samo mnożeń lub dzieleń. Koszt całej eliminacji wynosi więc 
252i (n — k)? m (2/3)n* operacji ary uaciyczny ch: Ponieważ układ 
trójkątny potrafimy rozwiązać kosztem n? (ASA koszt rozwiązania 
układu AZ = b jest również rzędu (2/3)n$ 


Opisany powyżej proces eliminacyjny załamie się gdy dla pewnego 
k element TSA 9 będzie zerem. W takim przypadku prosta modyfikacja 
polegająca na przestawieniu w k-tego układu z wierszem sķą-tym 
(k+1 < sk < n) dla którego aa + 7 0, pozwala kontynuować elimi- 
nację. Taki niezerowy element oczywiście istnieje, bo inaczej macierz 


AGD, a tym samym i macierz A, byłaby osobliwa. 
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Najczęściej wybieramy sy tak, aby 


co wymaga dodatkowo rzędu n*/2 porównań liczb rzeczywistych i żad- 
nych operacji arytmetycznych. Taki wybór elementu niezerowego pro- 
wadzi do eliminacji Gaussa z wyborem elementu głównego w kolumnie. 
Wtedy współczynniki ly = „rst (AR p. » są wszystkie co do modułu 
nie większe od jedności co, jak się przekonamy później, ma niebaga- 
telne znaczenie dla własności numerycznych całego algorytmu. 


3.3 Rozkład trójkątny macierzy 


Podamy teraz ważne twierdzenia wiążące eliminację Gaussa z rozkła- 
dem trójkątno-trójkątnym macierzy A. 


Twierdzenie 3.1 Eliminacja Gaussa beż przestawień wierszy (o ile 
jest wykonalna) jest równoważna rozkładowi macierzy A na iloczyn ma- 
cierzy trójkątnej dolnej L = (l; ;) (z jedynkami na przekątnej) i trójkąt- 
nej górnej R = (rij). Dokładniej, 


A = L-R, 
gdzie 
a, Ja; E 
lig 1 i= J, 
0 i<; 
oraz 


Dowód Zauważmy, że rugowanie elementów pod główną przekątną 
w k-tej kolumnie (k = 1,2,...,n — 1) odpowiada mnożeniu układu 
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AG-DE = -D z lewej strony przez macierz 


Ly = 1 
zli. 1 


= 1 


Ln—-1Ln-2 PSZ LaL} A =2 R 


AE "Li sa LiL 4) R. 


Łatwo sprawdzić, że macierz Lg" różni się od Ly tylko tym, że ele- 
menty —/,4 dla k +1 < i <n zmieniają się na l; x. Ponadto, mnożenie 
macierzy Lę" przez siebie w podanej kolejności odpowiada ” przepisy- 
waniu” elementów l; do macierzy wynikowej tak, że ostatecznie mamy 
PRZOD ES ANRE 


Zobaczmy teraz jak wybór elementu głównego w kolumnie i prze- 
stawianie odpowiednich wierszy macierzy układu wpływa na rozkład 
macierzy A z Twierdzenia 3.1. 

Przestawienie wiersza k-tego z s-tym w danej macierzy jest formal- 
nie równoważne pomnożeniu tej macierzy z lewej strony przez macierz 
permutacji elementarnej Py s. Jest to macierz, która powstaje z macie- 
rzy jednostkowej przez przestawienie jej k-tego wiersza z s-tym (albo 
k-tej kolumny z s-tą). Eliminacja Gaussa z wyborem elementu głów- 
nego w kolumnie odpowiada więc mnożeniu macierzy A z lewej strony 
kolejno przez macierze permutacji elementarnej Pķ,s„ i macierze *elimi- 
nacji” Ly tak, że 


EEE r E E ZOE | S E PAR: ewa LąP sa Lı Pi s, A = R. 
Ponieważ Pf, = I (identyczność), mamy 


PalhiPin A 6 Palai he Poa haA = LBP rA; 
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gdzie LU różni się od Lı tym, że zostały przestawione elementy (2, 1) 
i (52, 1). Następnie, 
Po sa LąLY” (Posa P11 A) 
(P3,ss Lo Ps...) (Ps, LI? Posa )(Po,ss Pasa Pis A) 
LL (Poea PrePiaiA), 


gdzie L? różni się od LV tym, że zostały przestawione elementy (3,1) 
i (s3,1), a LÊ? różni się od L+ przestawieniem elementów (3,2) i (83,2). 
Postępując tak dalej otrzymujemy 


R ENEA AR: bo PF, Lı PL; A 

DM "ole 1 | (Gole Mag) "HO 
gdzie DES” różni się od Ly jedynie pewną permutacją elementów w 
k-tej kolumnie pod główną przekątną, 

n—l1 
LÉ ) = L n-l,sn-1 ``’ Pk+1,sr41 LkPk+1,sr41 17 BR: ET. 

Postępując dalej tak, jak w przypadku eliminacji bez przestawień wier- 
szy, otrzymujemy następujący wniosek. 


Wniosek 3.1 Eliminacja Gaussa z wyborem elementu głównego w ko- 
lumnie jest wykonalna i jest równoważna rozkładowi 


PAETA 


gdzie P = Paisa tt Posa Pis, Jest macierzą permutacji, a L i R są 
macierzami konstruowanymi tak jak w Twierdzeniu 3.1, ale dla macie- 
rzy PA zamiast A. Ponadto wszystkie elementy macierzy L są co do 
modułu nie większe od jedności. © 


Dodajmy, że znalezienie czynników rozkładu PA = LR dla danej 
macierzy A kosztuje tyle samo co eliminacja Gaussa z wyborem ele- 
mentu głównego, czyli (2/3)n operacji arytmetycznych i n*/2 porów- 
nań. 

W praktycznej realizacji rozkładu PA = LR oczywiście nie musimy 
przy przestawieniu wierszy dokonywać fizycznie przepisań elementów 
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danego wiersza do innego. Wystarczy mieć dodatkowy n-wymiarowy 
wektor permutacji p|| interpretowany w ten sposób, że pļi] jest wskaź- 
nikiem do i-tego wiersza macierzy A. Przestawienie wierszy możemy 
wtedy realizować po prostu przez zamianę wskaźników, czyli elemen- 
tów wektora p. (Zob. U. 3.5.) 


Twierdzenie 3.1 i Wniosek 3.1 mają nie tylko teoretyczne znaczenie. 
Rozkładu macierzy na czynniki trójkątne (co w praktyce sprowadza się 
do pamiętania kolejnych permutacji i mnożników 1, ;) opłaca się doko- 
nywać w przypadku, gdy zadanie polega na rozwiązaniu nie jednego, 
ale wielu układów równań z tą samą macierzą A i ze zmieniającym się 
wektorem prawej strony, 


AŻ = b, dla 1<s<k. 


Jeśli każdy z tych układów rozwiązujemy *od początku” to musimy 
wykonać (2/3)n*k operacji arytmetycznych. Koszt ten możemy znacz- 
nie zredukować, jeśli dokonamy wstępnego rozkładu PA = LR. Wtedy 
układ Aż = b, jest równoważny układowi LRi = Pb., a więc jego 
rozwiązanie sprowadza się do rozwiązania dwóch znanych już układów 
trójkątnych: 

Lig = Pb, 


Taki sposób postępowania wymaga tylko (2/3)n* + 2n%k = 2n? (n/3+k) 
operacji arytmetycznych, co np. przy k = n powoduje obniżenie rzędu 
kosztu z ní do nè. 


Uwagi i uzupełnienia 


U. 3.1 Zaproponowany algorytm dla układów Ri = č z macierzą trójkątną 
R ma prawie optymalną złożoność kombinatoryczną (n?) wśród algorytmów 
dokładnych. Zauważmy bowiem, że jeśli algorytm jest dokładny to wykorzy- 
stuje wszystkie dane r,; i cy, 1 < i,j < n (zob. Ćw.3.1). Ponieważ każda 
kolejna operacja arytmetyczna (+, —, *, /) może wykorzystywać co najwyżej 
dwie nowe dane, a wszystkich danych jest n(n + 1)/2 + n z n*/2, algorytm 
dokładny musi wykonywać co najmniej n*/4 działań. 
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U. 3.2 Znane są algorytmy dokładne, które rozwiązują układ równań z ma- 
cierzą pełną kosztem proporcjonalnym do n”, gdzie p < 3, a więc mniejszym 
niż w eliminacji Gaussa. Algorytmy te są jednak nieużyteczne ze względu na 
skomplikowaną konstrukcję, zwykle dużą stłą przy n”, a także ze względu na 
złe własności numeryczne. Oczywiście, ograniczeniem dolnym na wykładnik 
p w koszcie rozwiązania układu z macierzą pełną jest 2 (bo mamy rzędu n? 
danych), ale nie wiadomo, czy może on być osiągnięty przez jakiś algorytm. 


U. 3.3 Układy równań z tą samą macierzą, ale ze zmieniającą się prawą 
stroną równania powstają często przy rozwiązywaniu, np. równań różniczko- 
wych cząstkowych, gdzie prawa strona układu odpowiada zmieniającym się 
warunkom brzegowym. Z wieloma układami tego typu mamy również do czy- 
nienia gdy chcemy znaleźć macierz odwrotną do danej macierzy A € R”*". 
Rzeczywiście, kolejne kolumny macierzy A”" są rozwiązaniami układów 


AŤ = Èj, 1<j<n, 
gdzie €; oznacza j-ty wersor. 


U. 3.4 Obok wyboru elementu głównego w kolumnie dokonuje się również 
wyboru elementu gównego w całej macierzy. To znaczy, w k-tym kroku eli- 
minacyjnym wybiera się element OWSA k<sp,tk < n, taki, że 


= max lafy 


Skok k<i,j<n l 


a następnie przestawia się wiersze k-ty z sąķ-tym i kolumny k-tą z ty-tą w 
macierzy A(*-1), Zauważmy, że przestawienie kolumn macierzy odpowiada 
mnożeniu jej z prawej strony przez macierz permutacji Pk tą. Przyjmując 
P= Pit Pota ttt Pn=1t,_, 1 korzystając z Twierdzenia 3.1 i Wniosku 3.1 
mamy, że eliminacja Gaussa z wyborem elementu głównego w całej macierzy 
jest równoważna rozkładowi 


PAP = LR. 


U. 3.5 Podamy teraz jedną z możliwych implementacji algorytmu rozkładu 
macierzy na iloczyn macierzy trójkątnych, PA = LR. W poniższym pro- 
gramie współczynniki macierzy wyjściowej są pamiętane w tablicy a|,|. Po 
wykonaniu programu, informacja o permutacji wierszy zapamiętana będzie 
w wektorze p|]. W alp[i],j] będą elementy (i, j) macierzy L dla i > j, oraz 
macierzy R dla i < j. 
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{ inicjacja wskaźników } 


for i:=l to n do pli] = i; 
for j:=1 to n—1 do 
begin 


{ wybór elementu głównego } 
im := j; val := abs(aļpļj], j]); 
for i:=ņj+1 to n do 


begin 
v := abs(a[p[i], j]); 
if v > val then 
begin 
im := i; val =i 
end; 


{ zamiana wskaźników | 
s := plim]; plim] := pl]; plj] := 8; 
{ eliminacja } 


v := als, j]; 
for i::j+1 to n do 
begin 
im := pli]; 
l := aim, j]/v; 
alim, j| := l; 
for k::j+1 to n do 
aim, k] := alim, k] — l» als, k] 
end; 
end; 
end. 


U. 3.6 Przedstawiony algorytm eliminacji Gaussa rozwiązuje układy z ma- 
cierzami dowolnej postaci. Jeśli macierz A jest szczególnej postaci to czasem 
proste jego modyfikacje pozwalają znacznie zmniejszyć koszt uzyskania roz- 
wiązania. Na przykład, niech A będzie macierzą trójdiagonalną postaci 


a1 C1 
bo a & 
A= b3 a3 cz 


bn ûn 
Wtedy układ A7 = € można rozwiązać stosując tzw. algorytm “przegania- 


nia”: 
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dh ::aq; ff i=e; 
for i:=2 to n do 
begin 
RE 
di := ai — l* cCi—1; 
fi := ei — lx fi 
end; 
11 := fn; 
for i:=n—1 downto 1 do 
cj = fi — i*r] 


Okazuje się, że jeśli macierz A ma dominującą przekątną, tzn. 
lail > [b] + |a|,  1<i<n, (3.3) 


(by = 0 = cn) i przynajmniej dla jednego i mamy powyżej nierówność <, 
to algorytm przeganiania jest wykonalny bez przestawień wierszy. Ponadto 
wymaga on 7n operacji arytmetycznych, a więc jest prawie optymalny. 


U. 3.7 Innym ważnym przykładem macierzy szczególnej postaci są macierze 
symetryczne i dodatnio określone. Są to macierze spełniające A = A” oraz 


dTAŻ>0, YZÆ0. 


Dla takich macierzy można nieco zmniejszyć koszt kombinatoryczny i zużycie 
pamięci przeprowadzając eliminację tak, aby otrzymać rozkład 


A= La Del? 


zamiast PA = LR, przy czym L jest tu jak zwykle macierzą trójkątną 
dolną z jedynkami na przekątnej, a D jest macierzą diagonalną z dodatnimi 
elementami na diagonali. 

Rozkład taki przeprowadzamy mnożąc macierz A przez znane już ma- 
cierze eliminacji nie tylko z lewej, ale też i z prawej strony, bez przestawień 
wierszy. Bez zmniejszenia ogólności rozpatrzymy tylko pierwszy krok. W 
tym celu, zauważmy najpierw, że 11 = ET AE, > 0 (gdzie €; jest pierw- 
szym wersorem), a więc nie musimy przestawiać wierszy, bo element na 
diagonali jest niezerowy. W pierwszym kroku mnożymy macierz A z lewej 
strony przez odpowiednią macierz Lı, a potem z prawej przez DŁ. Kluczem 
do zrozumienia algorytmu jest uwaga, że efektem mnożenia macierzy Lı A 
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z prawej strony przez LT? jest wyzerowanie elementów pierwszego wiersza 
poza a11 i pozostawienie niezmienionych pozostałych elementów. Ponadto 
macierz AU) = Lı ALI jest symetryczna i dodatnio określona. Rzeczywiście, 


(ADY" = (ALT = (LFYTATLĘ = LAL, 
oraz dla 1 Æ 0 
a" AQg = ZT L ALT? = (Li2)TA(LI2) > 0, 
bo z Æ 0 implikuje LTZ ź 0. Stąd ag = ET AE > 0. Postępując tak dalej 
otrzymujemy 
| PEYSIE TOO Loli ALT LE -LI LI; = D, 


przy czym macierz D jest diagonalna i dodatnio określona, a więc wyrazy 
na diagonali są dodatnie. Oznaczając 

L= Ly'Lq' -- Li, 
dostajemy żądany rozkład. 

Zauważmy, że przy praktycznej realizacji rozkładu A = LDLT wystarczy 
modyfikować jedynie wyrazy pod i na głównej przekątnej macierzy wyjścio- 
wej, ponieważ, jak zauważyliśmy, wszystkie kolejne macierze AG) są syme- 
tryczne. Pozwala to zmniejszyć koszt kombinatoryczny o połowę do n*/3 
operacji arytmetycznych. Opisany sposób rozkładu macierzy A = A” > 0 
nosi nazwę metody Banachiewicza-Choleskiego. (Zob. Ćw. 3.9.) 


U. 3.8 Rozkład trójkątno-trójkątny macierzy może być zastosowany rów- 
nież do policzenia wyznacznika macierzy A. Jeśli bowiem PA = LR to 


detA = (—1)*r1 1722:*:Tnn, 


gdzie s jest liczbą przestawień wierszy w eliminacji. 


Cwiczenia 
Ćw. 3.1 Pokazać, że jeśli algorytm rozwiązywania układu równań z ma- 
cierzą pełną (lub trójkątną) nie wykorzystuje wszystkich współczynników 


macierzy i wektora, to istnieją macierz A i wektor b takie, że algorytm dla 
tych danych nie daje dokładnego rozwiązania w arytmetyce idealnej. 


36 ROZDZIAŁ 3. UKŁADY RÓWNAŃ LINIOWYCH 


Ćw. 3.2 Znaleźć macierz odwrotną do macierzy trójkątnej. 


Ćw. 3.3 Pokazać, że iloczyn elementarnych macierzy permutacji jest ma- 
cierzą, która w każdym wierszu i w każdej kolumnie ma dokładnie jedną 
jedynkę. I odwrotnie, jeśli pewna macierz ma w każdym wierszu i kolumnie 
dokładnie jedną jedynkę to jest ona iloczynem pewnej liczby elementarnych 
macierzy permutacji. 


Ćw. 3.4 Pokazać, że jeśli dla danej macierzy A i permutacji P istnieje co 
najwyżej jeden rozkład PA = LR na macierz trójkątną dolną L z jedynkami 


na przekątnej i na macierz trójkątną górną R. 


Ćw. 3.5 Pokazać, że w algorytmie eliminacji Gaussa (bez przestawień) ele- 
menty l; i rij obliczane są według wzorów 


SE 
lij = (aig — XO Lima  j=1,2...,i— 1, 
k=1 


i—l1 
Tij = Qij — KSU? J=Li+1l,...,n, 
k=l 
dla i = 1,2,...,n. Wywnioskować stąd rozkład A = LR z Twierdzenia 3.1. 
Ćw. 3.6 Pokazać, że przy spełnieniu warunku (3.3) algorytm *przegania- 
nia” z U. 3.6 jest wykonalny bez przestawień wierszy. Opracować algorytm 
rozwiązywania układu z macierzą trójdiagonalną nie spełniającą warunku 


(3.3). 


Ćw. 3.7 Pokazać, że jeśli macierz pełna A ma dominującą przekątną, tzn. 
n 
2]a;.i| > DDICZJE 1<iśn 
j=l 


to jest to macierz dodatnio określona, 17 AZ > 0 dla £ Æ 0. Ponadto, elimi- 
nacja Gaussa jest dla takich macierzy wykonalna bez przestawień wierszy. 


Ćw. 3.8 Opracować algorytm eliminacyjny dla rozwiązywania układów li- 
niowych z macierzą Hessembega, czyli macierzą dla której a,; = 0, o ile 
i>j+2. 
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Ćw. 3.9 Opracować program rozkładający macierz symetryczną i dodatnio 
określoną A na iloczyn A = LDL" według metody Banachiewicza-Choleskie- 
go z U. 3.7. Obliczenia przeprowadzać w tej samej macierzy A pozostawiając 
elementy na i nad diagonalą nie zmienione i otrzymując w wyniku macierz 
L pod główną diagonalą macierzy A, a macierz D jako dodatkowy wektor. 


CE 
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ROZDZIAŁ 3. UKŁADY RÓWNAŃ LINIOWYCH 


Rozdział 4 


Analiza błędów w eliminacji 
Gaussa 


Zastanowimy się teraz, jak przebiega w arytmetyce fi, realizacja algo- 
rytmu eliminacji Gaussa z wyborem elementu głównego w kolumnie. 


4.1 Układy z macierzą trójkątną 


Najpierw pokażemy, że algorytmy rozwiązywania układów trójkątnych 
zaprezentowane w Rozdziale 3.1 są numerycznie poprawne. Rozpatrzmy 
najpierw układ trójkątny górny, Rx = c. Wykorzystując fakt, że licze- 
nie iloczynu skalarnego jest numerycznie poprawne, mamy 


fa = LEEA a 


1-0, | = 
Tnn(l RE zaa) ) Ba 
i dla i = n — 1,n — 2,..., 1, 
A,(z;) 


(a +6) — Sia rill + AAEN + %:5))(1 + wi) 
ria(l + Qii)(l + Ó,) 1 

c — Ljip Tiafa (e7) 

gdzie c = c;(1+5;), rh, = rall t aa) (1H, rż, = rigla) (1+ 

Yii) czyli |Ę — cl  v|cz| i |rż, — risl <(n+1)vlr,;|. Stąd otrzymane w 


r 
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fi, rozwiązanie jest dokładnym rozwiązaniem dla danych R” = (rę,)i.j 
i © =(c;);, przy czym 


IR” = R] < (n+ 1)v||Rllo i IE -= al < v |[E||oc. 
Zauważmy, że macierz zaburzona R” jest nieosobliwa. Powyższy algo- 
rytm jest więc numerycznie poprawny w klasie nieosobliwych macierzy 
trójkątnych górnych. 

Zauważmy również, że jeśli wektor prawej strony C jest reprezento- 
wany dokładnie (5, = 0), to błąd pochodzący z odejmowania można 
przenieść na zaburzenia r,, tak, że œ = č. Fakt ten wykorzystamy 
później. 

Przeprowadzając podobną analizę dla układu Li = č z macierzą 
trójkątną dolną z jedynkami na przekątnej dostajemy, że otrzymane w 
fi, rozwiązanie jest dokładnym rozwiązaniem układu L'x1 = œ, gdzie 

Z” — Zo S (+ 1)v||E]joo i Ie" = äle < © ||E||oo. 


4.2 Poprawność rozkładu trójkątnego 


Zajmiemy się teraz jakością numeryczną rozkładu macierzy na iloczyn 
trójkątno-trójkątny. W poniższym twierdzeniu norma macierzy jest do- 
wolna, ale ustalona. 


Twierdzenie 4.1 Niech A € R” będzie macierzą nieosobliwą. Dla 
dostatecznie silnej arytmetyki (tzn. dla dostatecznie małego v) rozkład 
trójkątno-trójkątny macierzy za pomocą eliminacji Gaussa z wyborem 
elementu głównego w kolumnie jest wykonalny w arytmetyce fi,. Otrzy- 
mane w wyniku macierze permutacji P”, trójkątna górna R” i trójkątna 
dolna L” (z jedynkami na przekątnej i elementami co do modułu nie 
większymi od jedności) spełniają równość 


Polas Boa daf 


gdzie 
|F|| < Kln)r||All, 


a K(n) jest pewną stałą niezależną odv i A. 
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Dowód Załóżmy najpierw dla uproszczenia, że eliminacja jest w fi, 

wykonalna i że nie musimy dokonywać przestawień wierszy, tzn. w 
; ; ; 33 (ASB) 

k-tym kroku elementem głównym jest obliczona wartość ay, , 


Zobaczmy jak zmieniają się elementy a ) przy realizacji algorytmu 
Gaussa w arytmetyce fl. Dla uniknięcia a indeksów, będziemy 
oznaczać przez w wartość w obliczoną w fl„. Rozpatrzymy dwa przy- 
padki, i <jii>jĵ. 

Niech najpierw i < j. Niech a >? będzie błędem bezwzględnym wy- 
(s) (s—1) 


tworzonym przy obliczaniu 2 aj; mając obliczone aj; , lis i 
4 (Błąd ten oszacujemy później). Baat 0 = a: ) —lik E 
otrzymujemy 
~ (k ~(k—1 7 „(k-11 k 
dj = al — lwdiy” + ej (4.1) 
- GRA - aga?) + E — aE 4 eE 
o adi (-1) | _(6) 
= Qij + | bi sås ; + eż, ) 
Wobec tego, że ÓW = aij, Mamy a® = = Qij + E e | < vlaizl. 


ży R = fij, bo ł-ty wiersz jest modyfikowany tylko w krokach 
od 1-go do (i — 1)-szego. Stąd, podstawiając k = i — 1, otrzymujemy 


i— 1 i—1 
e sli 2 (s) a 
Tij = Qij = lij + Den T a 
s=0 s=1 


Ponadto ać 


albo 2 
Qi j + ĉCij = Tij + DR (4.2) 
s=1 
gdzie e;;j = 555 z 
Dla i > j mamy podobnie 
Ad) _ „(s—1) (s) 
dj = dj + E F > (-i i adzy + R . (4.3) 


Wobec tego, że 


z 870 
lij = =) (1+ 0:5), 


IJ 
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otrzymujemy 
sd tx j ~(j—1 
af, Z li jà Gzy: 42. EJ < v|af; N, 
a stąd, (4.3) i z faktu, że a > ) = Sia 
Ig 
aij + eij = J lists j» (4.4) 
gdzie e;;j ZD, a 
Z (4.2) i (4.4) wynika równość 
(A+E)= L.R (4.5) 


z macierzą E = (e; j). Ponadto elementy A j Są co do modułu nie większe 
od jedności, ponieważ Eeh < jag >| (zob. Ćw. 1.2). 

Aby zakończyć dowód, trzeba jeszcze oszacować błąd |e, „|. W tym 
celu zauważmy, że dla k > 1 


-(k k-1 7 -(k-1 
w = (a | — wół (1 + w1)) (1 + wz), 


, 


Wą,wWą < v. Stąd i z (4.1) otrzymujemy 


ij 
k -(k -(k-1 7  _(k—1 ~ (k ai, ~ (k—1 
= wą af) (k—1) 
= Nad + wiri 
oraz 


(k) V alk) z(6—1) 
igl < yolli + lär |). 
Niech teraz Go = max; ; |a|, 


Gr = max [dy |, 
1) 


oraz G = MaXo<k<n-1 Gz. Dla i < j mamy 


Sl E G sa +G.)) 


1-2 


V V ; 
= NEM S r; i- DG, 


s=0 


IA 


lei j| 
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adlai >j 


A 


j 
esl < KEP < — ECES Gs-1 + G.)) 
s=0 


V V i 
x CEDNA Ar WEME 


W obu przypadkach 
2v 


1 — 


lei] < 


„TG. (4.6) 
Łatwo zauważyć, że a < (la;;| + lars) (1 + 20) oraz 
ail < däi I+ läk D +v), 
co implikuje Gy < 2G, 4, a w konsekwencji 
G<ZFT max |as]. 
Ostatnia nierówność oraz (4.6) dają 
e;| Ś n2”v max laij, 

a ponieważ dla dowolnej macierzy B € R”*” mamy 

-|Bl-» < maxļaij| < Bl. 

n ij 


to 
Ello < K(n)v |lAllæ, 


gdzie K (n) jest na poziomie n?2”. Z równoważności norm w R”*” wy- 
nika, że podobna nierówność zachodzi dla dowolnej, innej niż || : [loo 
normy, ale być może z inną stałą K(n). 

Oczywiście, ewentualne permutacje wierszy nie zmieniają tych osza- 
cowań, a powodują jedynie przemnożenie lewej strony równości (4.5) 
przez pewną macierz permutacji. 
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Aby zakończyć dowód, trzeba jeszcze pokazać, że dla dostatecznie 
małych v rozkład jest wykonalny, co jest równoważne temu, że macierz 
(A+ E) jest nieosobliwa. Rzeczywiście, weźmy v > 0 spełniające 


v K(n) |All Alo < 1. (4.7) 


Jeśliby macierz (A + E) była osobliwa to dla pewnego niezerowego 
wektora x mielibyśmy (4+E)x = A(I+A 'E)Z = 0, a w konsekwencji 
| ATEZ] = |lžllæ i 


1 < |A7'El|oo < |A7'||eol|F]|oo < K(nrlAllllA lo; 


co przeczy (4.7). 


Zauważmy, że warunek (4.7) na wykonalność rozkładu w fi, zależy 
od macierzy A poprzez wielkość 


cond(A) = |A|] IAT" 


Wielkość tą nazywa się uwarunkowaniem macierzy. Z Twierdzenia 4.1 
oraz warunku (4.7) wynika ważny wniosek o numerycznej poprawności 
rozkładu w klasie macierzy o uwarunkowaniu wspólnie ograniczonym 
przez pewną stałą. 


Wniosek 4.1 Niech M > 0. W klasie macierzy nieosobliwych A ta- 
kich, że 

cond(A) = |A A| < M 
rozkład PA = LR jest numerycznie poprawny. Dokładniej, dla v < 
1/(MK(n)) rozkład w fi, jest wykonalny, a otrzymane macierze P”, 
L” i R” pochodzą z dokładnego rozkładu macierzy (A+ E), gdzie || E|| < 
K(n)v|| A|. 


4.3 Poprawność rozwiązywania układu 


Przypomnijmy, że algorytm eliminacji Gaussa rozwiązywania układów 
równań AZ =b polega na rozkładzie macierzy na iloczyn PA = LR, a 
następnie na rozwiązaniu dwóch układów trójkątnych (albo, równoważ- 
nie, na sprowadzeniu układu do postaci trójkątnej i rozwiązaniu go). 
Numeryczna poprawność każdego z tych etapów algorytmu daje nam 
również numeryczną poprawność całego algorytmu. 
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Twierdzenie 4.2 Niech M > 0. Algorytm eliminacji Gaussa z wybo- 
rem elementu głównego w kolumnie zastosowany do rozwiązania układu 
równań 


Aż = b 
jest numerycznie poprawny w klasie danych wektorów b € R” oraz ma- 
cierzy nieosobliwych A E€ R”* spełniających 
cond(A) < M. 

Dowód Dla dostatecznie silnej arytmetyki, rozkład macierzy A daje 
w fi, P(A + E) = L Ry, gdzie ||E,|| < Kıv|| A||. W drugim kroku, 
otrzymane w fl, rozwiązanie jı układu Liy = Pıb jest dokładnym roz- 
wiązaniem układu (Lı + E>jy = Pb + €), gdzie ||Ez|| < Kəv|| L |l 
oraz ||e|| < Kavl|bl|. Z kolei w ostatnim kroku algorytmu, rozwiązując 
układ Rı£ = y, dostajemy w fl, rozwiązanie z”, które jest dokład- 
nym rozwiązaniem układu (Ry + E3)i = ty, gdzie ||E3|| < Kav||Rıll. 
(Zauważmy, że wektor ğı nie jest zaburzony, bo jest on reprezentowany 
dokładnie, zob. Rozdział 4.1). Stąd otrzymane w fl, rozwiązanie 7” jest 
dokładnym rozwiązaniem układu 

AD W, 
gdzie bb =b+E aA"=A+E spełnia równanie 

P(A + E) = (Lı + E2) (Ry + E3). 


Aby zakończyć dowód, należy teraz oszacować ||Æ||. Ponieważ P,(A + 
E) = Lı Rı, mamy 

PE = P E + L Es + Es Ro + E E3. 
Stąd 


|-F7||oo Ello + Za|ool|£3!|oo + l| E2 P2lls0 + || F2l|oc||-FE3||oo 
vKillAllo + ||E1|loo7 K3||Fa||so + I|Rillov KC || ||oo 


= v(KiljAll + (K2 + Ks)||Zul|ocll Ral|-c)- 


IA TA 


Ponieważ elementy macierzy Lı są co do modułu nie większe od jedno- 
ści, mamy ||Z4 ||ac < n. Jak zauważyliśmy wcześniej, mamy też || Ry ||c < 


46 ROZDZIAŁ 4. ANALIZA BŁĘDÓW W ELIMINACJI GAUSSA 


n2”-"||A||os. Przyjmując K = Kı + (K2 + K3)n*277! ostatecznie otrzy- 
mujemy 
lEl < Kv |All, 


czyli numeryczną poprawność, ponieważ stała K nie zależy od v i A. 
LI 


4.4 Uwarunkowanie macierzy, a błąd w 
fi, 


Pokazaliśmy, że eliminacja Gaussa jest numerycznie poprawna w klasie 
macierzy A € R”? takich, że 


cond(A) < M, 


gdzie M < oo jest jakąkolwiek stałą niezależną od A. Okazuje się, że 
wielkość uwarunkowania macierzy, cond( A), ma też zasadniczy wpływ 
na błąd zadania rozwiązywania układu równań. Rzeczywiście, mamy 
bowiem następujące twierdzenie. (Poniżej norma wektorowa jest do- 
wolna, ale ustalona, a norma macierzowa jest przez nią indukowana, 


zob. U. 1.2.) 


Twierdzenie 4.3 Niech E i € będą zaburzeniami odpowiednio macie- 
rzy A i wektora b takimi, że 


IE] < Kr] i lel < vlil, 


Jeśli 
Kivcond(A) < 1 


to układ zaburzony (A+ E)Z = (b+€) ma jednoznaczne rozwiązanie Ż* 
spełniające 


cond(A) 
1 — Kyvcond(A) 


|Z = 7| < (Ki + Kz)v ITI. 
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Dowód Zauważmy najpierw, że jeśli F jest macierzą taką, że |F|| < 1 
to macierz (I — F) (gdzie I jest macierzą identycznościową) jest nie- 


osobliwa oraz i 


rpe. 
Rzeczywiście, gdyby (I — F`) była osobliwa to istniałby niezerowy wektor 
z taki, że (I — F)z = 0, co implikuje |Fz||/||£|| = 1 i w konsekwencji 
IF] > 1. Aby pokazać (4.8) zauważmy, że 


IF=F) |< (4.8) 


1 = El=IF=PT=FP"| 
Le me a 


=E=], 


IV 


skąd bezpośrednio wynika (4.8). 


Po podstawieniu F = —A "FE mamy teraz 
IFI < IAHE < KiwlAll A] < 1, 


co wobec równości A + E = A(I + A™!E) daje, że macierz (A + E) 
jest nieosobliwa i układ zaburzony ma jednoznaczne rozwiązanie %*. 
Przedstawmy to rozwiązanie w postaci Z* = £* + (ż* — x"). Rozpisując 
układ zaburzony i wykorzystując równość Ax* = b otrzymujemy, że 
(A+E)(Z*-1*')=€ — Ex", czyli 

ZG = (I+ AE) A8- ET"), 
a stąd 


I*-=*"]| < E+ ATEA e + E 


AII > E 
< Kzv||bl| + Kiv || A|| lz" 
1 — Kywj| A|| |A- ( ) 
IAIA II A 
< (Kı + Kə) v ||t*||, 
1- Kiv|| All A= 
co kończy dowód. o 


Podsumowując Twierdzenia 4.2 i 4.3 otrzymujemy wniosek, który 
jest końcowym wynikiem tego rozdziału. 
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Wniosek 4.2 Niech A będzie macierzą nieosobliwą. Dla dostatecznie 
silnej arytmetyki, 
v K(n) cond(A) < 1 


(gdzie K(n) jest pewną stłą niezależną od A iv), eliminacja Gaussa z 
wyborem elementu głównego w kolumnie, zastosowana do rozwiązania 
układu AŻ = b, jest w fi, wykonalna i daje wynik fi,(1*) spełniający 
nierówność 


IA) — £*| < K(n)v cona(A) ||" |. 


Uwagi i uzupełnienia 


U. 4.1 Jak zauważyliśmy w dowodzie Twierdzenia 4.1 stała kumulacji K(n) 
zależy zasadniczo od wzrostu maksymalnego elementu Gy, w macierzach A) 
powstających w kolejnych krokach eliminacji. Okazuje się, że uzyskane, pe- 
symistyczne oszacowanie G = maxy Gy < 277!G5 jest praktycznie nie spo- 
tykane, choć teoretycznie może być osiągnięte. Przykładem jest macierz 


1 0 0- 0 1 
—1 1 0 : 1 
aie 4 : 1 
0 
1 
= el Zbeńss=f M 


dla której Ge = 2*Gy, 0 < k < n — 1, oraz G=2""!G5. 


U. 4.2 Eliminacja Gaussa bez przestawień wierszy jest numerycznie po- 
prawna, gdy element maksymalny w kolejnych macierzach A) wzrasta w 
sposób niezależny od A. Na przykład, gdy Gy < mGy_jtoG<m*"!Goiw 
Twierdzeniu 4.1 mamy K(n) = 2n?m*"!. Ma to miejsce wtedy gdy A jest 
symetryczna i dodatnio określona, albo ma dominującą główną przekątną, 
zob. U. 4.3 oraz Ćw. 4.5 i 4.6. 


U. 4.3 Dla macierzy symetrycznych i dodatnio określonych, A = A” > 0, 
eliminacja Gaussa jest wykonalna bez przestawień wierszy, zob. U. 3.7. W 
klasie tych macierzy eliminacja bez przestawień wierszy jest też numerycz- 
nie poprawna, a więc w szczególności algorytm Banachiewicza-Choleskiego 
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jest numerycznie popeawny. Wykażemy to pokazując, że maksymalny ele- 
ment w kolejnych macierzach A) wzrasta co najwyżej dwukrotnie, tzn. 
Gk < 2Gk-1, 1 < k < n. W tym celu wykorzystamy fakt, że dla dowolnej 
symetrycznej i dodatnio określonwj macierzy B = (b,;) mamy b, > 0, oraz 
spełniona jest nierówność 


b; < b ibj j, Vi, j (4.9) 


(zob. Ćw. 4.1). Jak wykazaliśmy w U. 3.7, każda z macierzy A) jest syme- 
tryczna i dodatnio określona, Stąd 


k k— 
lagj| = lafy? — zaj, | 
(-) Ense E 
< l *|+ ylk] 
ae] 

(4-1) (k=1) (k-1) 
< Ge M 1) j Vaii Va Akk jj 
= td aj GD 

a 3 
(k—1) (k—1) (k— 1) (k—1) 
2/aj; ajj < Zmax(aj; aj; D 


a w konsekwencji Gy < 2Gy_11. 


Cwiczenia 


Ćw. 4.1 Wykazać, że macierz symetryczna 2 x 2 


a c 

c b 
jest dodatnio określona wtedy i tylko wtedy, gdy a > 0 i ab > œ. Wywnio- 
skować stąd nierównosc (4.9) dla macierzy symetrycznych i dodatnio określo- 


nych dowolnego wymiaru. Ponadto, największy co do modułu element takiej 
macierzy leży na głównej diagonali. 


Ćw. 4.2 Wykazać, że macierz A jest dodatnio określona wtedy i tylko wtedy 
gdy dla każdego i € R” wektory Ar i £ tworzą w R” kąt ostry. 
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Ćw. 4.3 Pokazać, że jeśli eliminację Gaussa z wyborem elementu głównego 
w kolumnie zastosujemy do macierzy trójdiagonalnej, to wzrost elementu 
maksymalnego macierzy nie będzie zależał od n. Dokładniej, 


(k) X 
max |a; | <2 max |as j]. 


Ćw. 4.4 Pokazać, że dla macierzy Hessenberga (aij = 0 dla i > j +2) 
eliminacja Gaussa z wyborem elementu głównego w kolumnie daje 


max jaf)| < (k + 1) max |a; j|. 
i j,k wJ 


Ćw. 4.5 Pokazać numeryczną poprawność algorytmu przeganiania z U. 3.6. 


Ćw. 4.6 Wykazać numeryczną poprawność eliminacji Gaussa bez przesta- 
wień wierszy dla macierzy z dominującą przekątną, tzn. gdy 


n 


2|a;,;| > DDICZJE l<i<n. 
j=1 
z A (b) E D 
Wskazówka. Zauważyć, że max;,j |a; | < 2 max; j |a; | 
Ćw. 4.7 Jeśli ; 
(A+ E)? = b, (4.10) 


gdzie |F||, < Kv|AJ|p, to oczywiście dla residuum ř = b— AZ mamy 
[Flp < KrA|pliz||p- (4.11) 


Pokazać, że dla p = 1,2,00 zachodzi też twierdzenie odwrotne, tzn. jeśli 
spełniony jest warunek (4.11) to istnieje macierz pozornych zaburzeń E taka, 
że ||F||p < Kv||A||p oraz spełniona jest równość (4.10). 

Wskazówka. Rozpatrzyć E = r(sgn ż;)f-;e,/||Z]i da p = 1, E=rz"/||z||3 
dla p = 2, oraz E = r(sgn zę)€k /||Z||oc dla p = oo, gdzie k jest indeksem dla 
którego |z| = ||Z||oo. 


Rozdział 5 


Zadanie wygładzania 
liniowego 


W tym rozdziale zajmiemy się zadaniem wygładzania liniowego, na- 
zywanym też liniowym zadaniem najmniejszych kwadratów. Jest ono 
uogólnieniem zadania rozwiązywania kwadratowych układów równań 
liniowych do przypadku, gdy układ jest nadokreślony. 


5.1 Układ normalny 


Niech A będzie daną macierzą o m wierszach i n kolumnach, A € R”*”, 
taką, że 
m > n = rank(A), 


albo równoważnie, taką że jej wektory kolumny są liniowo niezależne. 
Niech także dany będzie wektor b € R”. Jasne jest, że wtedy układ 
równań AZ = b nie zawsze ma rozwiązanie - mówimy, że układ jest 
nadokreślony. 

Zadanie wygładzania liniowego polega na znalezieniu wektora £* € 
R”, który minimalizuje wektor residualny r = b— AT w normie drugiej, 
tzn. A 3 

|Ë - A£|» = min [6 — Az. 
xeR" 


Przykład 5.1 Przypuśćmy, że dla pewnej funkcji f : [a,b] + R ob- 
serwujemy jej wartości f; (dokładne lub zaburzone) w punktach t;, 
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1 < i < m. Funkcję tą chcielibyśmy przybliżyć inną funkcją w nale- 
żącą do pewnej n wymiarowej przestrzeni liniowej W, np. przestrzeni 
wielomianów stopnia mniejszego niż n. Jakość przybliżenia mierzymy 
wielkością 


(di - wlt)? (5.1) 


Wybierając pewną bazę (w;);_, w W i rozwijając w w tej bazie, w(t) = 
2j=1 Cjw;(t), sprowadzamy problem do minimalizacji 


m 


2 
S(r- Fee) 
i=1 j=1 
względem cj, a więc do zadania wygładzania liniowego. Rzeczywiście, 
kładąc A= (aij) ER”? z Qij = w;(t;), b = (WEG LG = (cj )j=1, 
wielkość (5.1) jest równa ||b — AŻ||3. 


Lemat 5.1 Zadanie wygładzania liniowego ma jednoznaczne rożwiąża- 
nie £*, które spełnia układ równań 


ATAG = ATb. (5.2) 


Dowód Niech P C R” będzie obrazem A jako odwzorowania liniowego 
z R” w R”, 
P = { AT: ZER”}. 

Ponieważ kolumny macierzy A są liniowo niezależne, tworzą one bazę 
w P. Stąd dim(P) = n i odwzorowanie A : R” — P jest różnowar- 
tościowe. Ponadto przestrzeń R™ z normą drugą jest przestrzenią uni- 
tarną. Residuum jest więc minimalizowane dla wektora £* € R” ta- 
kiego, że AT* jest rzutem prostopadłym wektora b na podprzestrzeń P. 
(Przypomnijmy, że skończony wymiar P zapewnia, że taki rzut istnieje 
i jest wyznaczony jednoznacznie.) Równoważnie można powiedzieć, że 
residuum b — A7* jest prostopadłe do P, 


(AZ)T(b— A7") = 0, YZER”, 


albo e 
rT (ATb — ATAT) = 0. YZeER. 
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Otrzymaliśmy więc, że wektor ATB — AT AT* jest prostopadły w R” do 
każdego innego wektora. Ponieważ jedynym wektorem o tej własności 
jest wektor zerowy, to AT AT* = ATb. 


Zauważmy, że jeśli macierz A jest kwadratowa, m = n, to rozwią- 
zaniem jest 7* = Ab i residuum jest zerem. Zadanie wygładzania li- 
niowego jest więc uogólnieniem rozwiązywania kwadratowych układów 
równań liniowych. 

Równanie (5.2) nazywa się układem normalnym. Może ono nam 
sugerować sposób rozwiązania zadania wygładzania liniowego. Wystar- 
czy bowiem pomnożyć macierz A” przez A i rozwiązać układ normalny. 
Zauważmy ponadto, że macierz AT A jest symetryczna i dodatnio okre- 
ślona, bo (A7A)7 = ATA i dla Z # 0 mamy £7(A7A)7 = (AZ)7”(AŻ) = 
|Az||3 > 0, przy czym ostatnia nierówność wynika z faktu, że kolumny 
macierzy A są liniowo niezależne i dlatego A£ +Æ 0. Przy mnożeniu A” 
przez A wystarczy więc obliczyć tylko elementy na głównej przekątnej 
i pod nią, a do rozwiązania równania z macierzą A” A można zastoso- 
wać algorytm Banachiewicza-Choleskiego opisany w U. 3.7. Jak łatwo 
się przekonać, koszt takiego algorytmu wynosi n*(k + n/3), przy czym 
dominuje koszt mnożenia obliczenia macierzy ATA. 

Ma on jednak pewne wady. Mnożenie macierzy powoduje w fi, po- 
wstanie “po drodze” dodatkowych błędów, które mogą nawet zmienić 
rząd macierzy. Na przykład, dla macierzy 


1 1 1 1 


mamy 


2 


€ 


1 
1 
1 1+e? 1 
1 1 1+ 


Jeśli £? < v to fl,(1+£°) = 1, co implikuje rank(fl, (AT A)) = 1, podczs 
adya A 4 
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Poniżej przedstawimy inną metodę rozwiązywania zadania wygła- 
dzania liniowego, która oparta jest na specjalnych przekształceniach 
zwanych odbiciami Hauseholdera. 


5.2 Odbicia Hauseholdera 
Dla danego wektora w € R” o normie |w|- = vw7w = 1, odbicie 
(macierz) Hauseholdera zdefiniowane jest jako 


Zauważmy, że 


a ponieważ (W”r)w = (Z, W)aw jest rzutem prostopadłym 7 na kieru- 


nek wektora w ((:,:)2 oznacza iloczyn skalarny), to Hg jest odbiciem 
lustrzanym wektora 7 względem hiperpłaszczyzny (wymiaru m — 1) 
prostopadłej do w. 
Odbicia Hauseholdera są przekształceniami nieosobliwymi spełnia- 
jącymi 
H` = H = H”. 
Rzeczywiście, ponieważ w ma normę jednostkową, mamy 


H? = (I — 200") = I — 4d” + 4u(wujw? = I, 


Oraz 
H" = (I -200 = I- Uà) w = I. 


W szczególności H jest więc przekształceniem ortogonalnym, H7} = 
HT, czyli nie zmienia długości wektora, 


|Hz|- = y(H2)"(H2) = JaT(HTH)A = VITZ = |a|». 


Odbicia Hauseholdera zastosujemy do przeprowadzenia danego wek- 
tora a Æ 0 na kierunek innego niezerowego wektora, powiedzmy €, tzn. 
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Załóżmy dla uproszczenia, że ||e||2 = 1. Aby wyznaczyć H zauważmy, 
że 


a ponieważ a = +||a||> i ||W||2 = 1 to 


GF Ilall 
la 7 |d||>e]|2 


W szczególności, jeśli € = €, jest pierwszym wersorem to powyższe 
wzory dają 


——T 
Fa 2 Gaj 
yY 
gdzie 
_ ] aFläl i=1, 
U; = > 
ai Ż<Xi<m, 
oraz 


i 1 p m 
y = zla = 5((a F lal) + 37 aż) 
i=2 
1 i 2 >||2 2 
= s(Ż2a + [alk = allele) = lal = alal 
i=1 


Otrzymaliśmy dwa odbicia Hauseholdera przekształcające dany wektor 
d na kierunek pierwszego wersora, w zależności od wybranego znaku 
przy ||a||z. Ustalimy ten znak na plus gdy a, > 0, oraz na minus gdy 
aq < 0, co pozwoli na obliczenie w i y z małym błędem względem w 
fi. Wtedy bowiem mamy 


wwe a1 + |a||> ai > 0, 
aj — |a||> Qı < 0, 


oraz y = ||a||3 + |a| ||d||2, czyli zawsze dodajemy liczby tych samych 
znaków. Ponadto pierwsza współrzędna wektora Ha jest równa — ||a||> 
dla a, > 0, oraz +-||a|| dla a, < 0. 
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5.3 Algorytm dla zadania wygładzania 


Odbić Hauseholdera można użyć do rozkładu macierzy A € R””" na 
iloczyn ortogonalno-trójkątny. 

Niech A = (dy, d2, . . . , En), gdzie d; są wektorami-kolumnami macie- 
rzy A. Wybierzmy pierwsze odbicie Hauseholdera Hy = Im — wt? /y1 
tak, aby przekształcało pierwszy wektor-kolumnę macierzy A na kieru- 
nek e,. Efektem pomnożenia macierzy A z lewej strony przez Hı będzie 
wtedy macierz 


AO = (a®,...,8®) = (Hyfh,..., Hd), 


w której pierwsza kolumna af) ma niezerową tylko pierwszą współ- 


rzędną. W następnym kroku wybieramy drugie przekształcenie Hause- 
holdera Hə = Im—1 — DUL /%2 wymiaru m — 1 tak, aby przeprowadzało 
wektor (aj; Ji. na kierunek pierwszego wersora w R”*"'. Rozszerzając 
üə E€ R”! do wektora ile € R” przez dodanie zera jako pierwszej 
współrzędnej, tlo = (0, 02)”, otrzymujemy przekształcenie (macierz) 
Hauseholdera Ho = Im — Ttk /Y2 w R” postaci 


n 
Beee h 
0 Ho 


Pomnożenie macierzy AU) z lewej strony przez Ha spowoduje teraz 
wyzerowanie drugiej kolumny macierzy pod elementem Ga; przy czym 
pierwszy wiersz i pierwsza kolumna pozostaną niezmienione. Postępując 
tak dalej n razy (albo n — 1 razy gdy m = n) otrzymujemy 


H,„Hn-1: HHA = R, 

gdzie R € R”** jest uogólnioną macierzą trójkątną górną, tzn. r;j; = 0 
dla ¿ > j. Stąd, podstawiając Q = H,H>---H,, dostajemy rozkład 
macierzy na iloczyn ortogonalno-trójkątny 

A=Q:R. (5.3) 
Rzeczywiście, macierz Q € R™*™ jest ortogonalna, bo 

Q! = (HHo:--H„) ' =H.!---H;"'H;" 
= H... HEH? = (HiH; Ha) = Q". 
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Dyspunując rozkładem (5.3) zadanie wygładzania liniowego można 
rozwiązać następująco. Ponieważ mnożenie przez macierz ortogonalną 
nie zmienia normy drugiej wektora, mamy 


IF = |b- Az] = |ë- QR2|- 
IQ(Q7?- RZ)|» = |E— Rz», 


gdzie c = Q7Tb = = -- HąHqb. Rozbijając wektor € na Z = (ĉr, rr)”, 
gdzie ĉr € R” i ĉr > R”-”, oraz macierz R na 


Rej k 


gdzie Ry € R”*%” jest macierzą trójkątną górną, a 0 jest macierzą zerową 
wymiaru (m — n) X n, otrzymujemy 


Irl = IE- Faź||z + Err 


Rozwiązanie £* zadania wygładzania jest więc rozwiązaniem układu 
liniowego trójkątnego, 


—: 


AR R;'ćh, 
oraz |i*|> = ||5— Aż*| = ||ćzr|>. 


Zastanówmy się nad praktyczną realizacją tego algorytmu. Każde z 
kolejnych przekształceń Hauseholdera Hy wyznaczamy przez obliczenie 
yk oraz współrzędnych wektora tę. Wektor ten pe tylko m-k+1 
współrzędnych niezerowych, a ponadto u; = = A ) dla k+1 SEM: 
Dzięki takiej reprezentacji Hy, mnożenia Hy możemy dla dowolnego 
x realizować według wzoru 


(HyZT); = Ti — SUki, 


gdzie s = uj £/4%. 

Uwzględnizjąc obecność zerowych elementów w u, przejście od ma- 
cierzy AFM do AW) kosztuje rzędu 4(m—k+1)(n— k) operacji arytme- 
tycznych i obliczenie jednego pierwiastka kwadratowego. Cały rozkład 
A = QR kosztuje więc rzędu (dla dużych m i n) 


),4(m— k+1)(n=k) m Sn? + 2n?(m — n) = 2n'(m — n/3) 
k=1 
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operacji arytmetycznych i n pierwiastków kwadratowych. Zauważmy, 
że w przypadku m = n, a więc dla kwadratowego układu równań, koszt 
ten wynosi (4/3)n* i jest dwa razy większy od kosztu eliminacji Gaussa. 


Uwagi i uzupełnienia 


U.5.1 Pokazaliśmy, że rozwiązaniem zadania wygładzania liniowego jest 
wektor £* = A*b, gdzie 


At = (474)""A7. 


Macierz A* € R”*™ nazywa się macierzą pseudoodwrotną do A € R”*", 
rank(A) = n. Dla nieosobliwych macierzy kwadratowych mamy oczywiście 
AT=A"!, ponieważ wtedy f* = A |. 


U. 5.2 Pokażemy, że każdą macierz A”*", rank(A) = n, można rozłożyć 
na iloczyn 


A=UxzV7, (5.4) 


gdzie U € R™®%™ i V e R”*” są macierzami ortogonalnymi, a X € R™%” 
jest macierzą diagonalną (tzn. o;i; = 0 dla i £ j) o dodatnich wyrazach øg; ;. 

Ponieważ macierz ATA jest symetryczna i dodatnio określona, znane 
twierdzenie z algebry liniowej mówi, że istnieje w R” baza ortonormalna 
( 4 wektorów własnych tej macierzy, a odpowiadające im wartości własne 
są rzeczywiste i dodatnie, tzn. 


oraz (ATAJE; = N;ćj, gdzie 
A2 A22- d> An 70. 
Zauważmy, że wektory 1; = À; 1/ (AE) są ortonormalne w R”, bowiem 


(aoti = (AA) AE, AE) 
(AA) PUATA), &)2 


= (AA TANE, E) = i > siak, 
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Wektory tę, 1 < i < n, można uzupełnić m — n dodatkowymi wektorami 
tak, aby cały układ (m): był bazą ortonormalną w R”. Zdefiniujmy teraz 
macierze ortogonalne U o kolumnach h, 1 < i < m, oraz V o kolumnach 
E, 1 < j < n. Bezpośredni rachunek pokazuje, że macierz UT AV jest dia- 
gonalna z wyrazami na przekątnej ,/A;. Stąd A = UXVT z Ojj = yA 
T1EJŚm 

Liczby „À; nazywa się wartościami szczególnymi macierzy A, a rozkład 
(5.4) rozkładem macierzy według wartości szczególnych. 

Zauważmy jeszcze, że jeśli A = UZV”T to At = V>TU? oraz Xt € 
R”*™ jest macierzą diagonalną z wyrazami na diagonali równymi 1//N;. 
Rozwiązanie zadania wygładzania liniowego można więc zapisać jako 


i = vVvxntUTb. 


U.5.3 Przedstawimy teraz jedną z możliwych implementacji algorytmu 
Hauseholdera rozkładu macierzy na iloczyn A = QR. Po wykonaniu poniż- 
szego programu wyrazy rij macierzy R dla 1 < 4 < j zostaną zapamiętane 
na miejscach ali, j], współrzędne uj; wektora ii; dla j +1 < 4 £ m na miej- 
scach ali, k], a współrzędna uj; i liczba y; odpowiednio na ulj] i gam[j], 
LSJŚmM 


for k:= 1 to n do 


begin 
4 obliczanie kolejnego odbicia Hauseholdera | 
norm2 := 0.0; 
for l:k to m do 
begin 
aa := ajl, k]; 
norm2 := norm2 + aa * aa 
end; 
norm := sqrt(norm2); 
aa := ajk, kl]; 
if (aa > 0.0) then 
begin 
uu := aa + norm; 
akk := —norm 
end else 
begin 
uu := aa — norm; 


akk := norm 
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end; 
gamma := norm2 + abs(aa) * norm; 
ujk] := uu; 
gamjk]| := gamma; 
{ modyfikacja kolumn macierzy } 
alk, k] := akk; 
for j:k+1 to n do 
begin 
s := uu x alk, j]; 


for l:=k+1 to m do 
s := s + all.k| x afl, jl; 

s := s/gamma; 

a|k,j] := alk,j] — s * uu; 

for l:k+1 to m do 
all, j| := ajl, j| — s*all,k] 

end; 
end. 


U. 5.4 Można pokazać, że przedstawiony algorytm Hauseholdera rozkładu 
macierzy na iloczyn ortogonalno-trójkątny jest numerycznie poprawny. To 
znaczy, otrzymane w fl, macierz trójkątna górna R” i ortogonalna Q” (re- 
prezentowana przez n wektorów ty i liczb %,) spełniają (A + E) = Q” FY, 
gdzie | E|| < K(m,n)v|All 


U.5.5 Zadanie wygładzania liniowego można również rozwiązać stosując 
innego rodzaju rozkład ortogonalno-trójkątny macierzy A € R™*”; miano- 
wicie przez zastosowanie ortogonalizacji Grama-Schmidta do wektorów ko- 
lumn da; macierzy A. W wyniku otrzymujemy ciąg wektorów q,, 1< j <n, 
tworzący układ ortonormalny w R”, oraz spełniający 


spanfqi,q>,...,qjł = spanfdy,do,...,d;), 1<j<n. (5.5) 


Zauważmy, że jeśli z wektorów q; stworzymy macierz Q = (q1,fb,..-;dn) € 
R”, to (5.5) implikują istnienie kwadratowej macierzy trójkątnej górnej 
Re R”? takiej, że 

A=Q-R. 


(W odróżnieniu od efektu działania algorytmu Hauseholdera, macierz Q nie 
jest tu kwadratowa, ale za to kwadratowa jest macierz R.) 
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Wektory ortonormalne q; oraz współczynniki r, ; macierzy R można wy- 
znaczyć z układu równań 


j 
dj = X rs jls, 1<j<n. 
s=l 


Mnożąc j-te równanie skalarnie przez q,, 1 < i < j, otrzymujemy fij = 
(dj, qi)2 = qq dj, a stąd wzory rekurencyjne 


for j := 1 to n do 


begin 
for i:=1 to j—1 do rij := (Ẹ,@j)2; 
P J=L uS, 
Dj = aj — 251 Ts jls; 
rjj = |Bzl|2; 
Ü = Pi/Tj 
end. 


Dysponując rozkładem A = QR, rozwiązanie x* zadania wygładzania 
wyznaczamy z równania Ri = QTb. Rzeczywiście, układ (G)j=1 można for- 
malnie uzupełnić do układu ortonormalnego w R” dodając do niego pewne 
wektory q; dla n+1 < j < m. Tworząc macierz ortogonalną Q = (q,..., dm) 
i rozszerzając macierz R € R”? do macierzy R € R™*” przez dodanie do 
niej (m—n) zerowych wierszy, otrzymujemy A = QR, a więc rozkład taki jak 
w algorytmie Hauseholdera. Postępując dalej tak, jak w analizie algorytmu 
Hauseholdera, otrzymujemy żądany wynik. 


U. 5.6 Okazuje się, że algorytm ortogonalizacyjny Grama-Schmidta z U. 
5.5 ma niedobre własności numeryczne, gdy wektory-kolumny macierzy wyj- 
ściowej A są “słabo liniowo niezależne”; tzn. otrzymane w fi, wektory T 
mogą być wtedy “słabo” ortogonalne. W takim wypadku należy stosować 
podwójną ortogonalizację, najpierw do wektorów d;, a potem do TE 1< 


j < n. (Zob. Ćw. 5.7.) 


Cwiczenia 


Ćw. 5.1 Pokazać, że dla macierzy pseudoodwrotnej A* do danej macierzy 
A € R”*", rank(A) = n, macierz ATA = I, jest identycznością w R"*", 


natomiast 
In 0 
+— n mxm 
AA = ( 0 0 | ER ; 
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czyli AA? jest rzutem prostopadłym na podprzestrzeń rozpiętą na pierw- 
szych n wersorach w R”. 


Ćw. 5.2 Niech A € R”*7. Uzasadnić, że jądro macierzy A”, 
ker(A”) = {gE R”: ATg=0), 

jest podprzestrzenią prostopadłą do obrazu 4, 
im(4) = { ATER”: żeR”). 


Ćw. 5.3 Pokazać, że macierze AT A i AAT mają takie same niezerowe war- 
tości własne, a podprzestrzenie własne im odpowiadające mają ten sam wy- 
miar. 


Ćw. 5.4 Uzasadnić, że dana macierz kwadratowa Q € R™*™ jest ortogo- 
nalna wtedy i tylko wtedy gdy jej kolumny (wiersze) tworzą bazę ortonor- 
malną w R”. 


Ćw. 5.5 Niech ii będzie niezerowym wektorem w R™ oraz y = [u||3/2. 
Uzasadnić, że algorytm obliczania 


Hi = 7 — si, s = (W z)/y, 


według powyższego wzoru, jest numerycznie poprawny ze względu na dany 
wektor i € R”. 


Ćw. 5.6 Policzyć złożoność algorytmu ortogonalizacyjnego opisanego w U. 
5.5 rozwiązania zadania wygładzania liniowego. 


Ćw. 5.7 Zastosujmy ortogonalizację Grama-Schmidta do dwóch wektorów 
liniowo niezależnych G i b o normach dl» = 1 = ||b||2. Załóżmy dla uprosz- 
czenia, że w obliczeniach jedynie iloczyn skalarny tych wektorów s = (a, DE 
liczy się z błędem, fi,(s) = s(1 + £), gdzie |e| jest dodatnie i na poziomie 
v. Pokazać, że wtedy dla otrzymanej w wyniku ortogonalizacji unormowanej 
pary wektorów a i č mamy 


TES 
1 — s21- 


(a, €)> = 


Wywnioskować stąd, że gdy ai b są “prawie liniowo zależne”, to aicC są 
dalekie od ortogonalnych. 
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Ćw. 5.8 Zaprogramować algorytm Grama-Schmidta z U. 5.5 (z podwójną 
ortogonalizacją) rozwiązujący zadanie wygładzania. 


Ćw. 5.9 Zaprogramować algorytm obliczający macierz odwrotną do danej 
nieosobliwej macierzy A € R”*? wykorzystujący rozkład Hauseholdera A = 


QR. 
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Rozdział 6 


Interpolacja wielomianowa 


Dotychczas rozpatrywaliśmy zadania, w których danymi są ciągi liczb 
rzeczywistych. Teraz zajmiemy się zadaniami, w których danymi są 
funkcje o wartościach rzeczywistych. Pierwszym z nich jest zadanie in- 
terpolacji wielomianowej. 


6.1 Sformułowanie zadania interpolacji 


Niech D C R i niech F będzie pewnym zbiorem funkcji f : D > R. 
Niech £o, 11,... , £n będzie ustalonym zbiorem parami różnych punktów 
z D, zwanych później węzłami. 

Powiemy, że wielomian w interpoluje funkcję f € F w węzłach z;,, 


gdy 
w(u;) = f(z;  0<j<n. 


Oznaczmy przez II, przestrzeń liniową wielomianów stopnia co naj- 
wyżej n o współczynnikach rzeczywistych, 


In, = { w(x) = ans” +a it" TE" kati GGERÓLJ<NE 
j 


Lemat 6.1 Dla dowolnej funkcji f : D > R istnieje dokładnie jeden 
wielomian w; € Il, interpolujący f w węzłach tj, 0<j<n. 


Dowód Wybierzmy w II, dowolną bazę wielomianów g,, 0< j <n, 
II = span{ Po, Pls- - -Pn F 
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Wtedy każdy wielomian z II, można jednoznacznie przedstawić w po- 
staci rozwinięcia względem wybranej bazy. Warunkiem koniecznym i 
dostatecznym na to, aby wielomian w;(:) = X5-6 G;9;(:) interpolował 
f jest spełnienie układu n + 1 równań liniowych 


z n+ 1 niewiadomymi c,, który w postaci macierzowej wygląda nastę- 
pująco: 


po(to) Pi(to) * Pn(zo) Co f (zo) 
polxı) pılzı) * Pn(11) C1 E „> l (6.1) 
Po(Tn) Pi(Tn) > Pn(Tn) Cn f(n) 


Aby wykazać, że układ ten ma jednoznaczne rozwiązanie wystarczy, aby 
wektor zerowy był jedynym rozwiązaniem układu jednorodnego. Rze- 
czywiście, układ jednorodny odpowiada interpolacji danych zerowych, 
f(c;) = 0, Vi. Istnienie niezerowego rozwiązania byłoby więc równo- 
ważne istnieniu niezerowego wielomianu stopnia nie większego od n, 
który miałby n + 1 różnych zer z;, co jest niemożliwe. 


Zadanie znalezienia dla danej funkcji f jej wielomianu interpolacyj- 
nego stopnia co najwyżej n jest więc dobrze zdefiniowane, tzn. rozwią- 
zanie istnieje i jest wyznaczone jednoznacznie. Zauważmy, że wielomian 
interpolacyjny wy jako taki nie może być wynikiem obliczeń w naszym 
modelu obliczeniowym, możemy natomiast wyznaczyć jego współczyn- 
niki cj w wybranej bazie. 


Definicja 6.1 Niech (,p;);_, będzie bazą w przestrzeni Il„ wielomianów 
stopnia co najwyżej n. Zadanie (obliczeniowe) interpolacji wielomino- 
wej polega na obliczeniu dla danej funkcji f współczynników c; takich, 
że wielomian 


wl) = Keil) (6.2) 


interpoluje f w punktach z,,0<j<n. 
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6.2 Wybór bazy wielomianowej 


Można powiedzieć, że trudność zadania interpolacji jest zdetermino- 
wana przez trudność rozwiązania układu (6.1), ta zaś zależy w istotny 
sposób od wybranej bazy (,;)7-9. W naturalny sposób powstaje więc 
problem wyboru “wygodnej” bazy w II,. Rozpatrzymy trzy bazy: La- 
grange'a, potęgową i Newtona. 


BAZA LAGRANGE'A (KANONICZNA) 
Zdefiniujmy dla 0 < j < n wielomiany 
(z — zo) (z = 1): +: (£ — 7-1 )(E — zaa) +: (E — 21) 


(£; — ©o)(2; — £1) + (Ej — Lj-1)(Tj — Lj+1) (2j — Zn) 
(6.3) 


(z) = 
Zauważmy, że każdy z l, jest stopnia dokładnie n oraz 
[o ižj, 
as l 1 i=j 


Stąd łatwo widać, że wielomiany te stanowią bazę w IIn, którą nazy- 
wamy bazą Lagrange’a. Macierz układu (6.1) jest identycznością i w 
konsekwencji c; = f(x;), Vj. Wielomian interpolacyjny dla funkcji f 
można więc w tym przypadku zapisać jako 


wę() = X fC). 
j=0 
Koszt kombinatoryczny rozwiązania zadania interpolacji jest przy tym 
zerowy. 
Przypuśćmy, że chcielibyśmy obliczyć wartość wielomianu interpo- 
lacyjnego w; w punkcie z różnym od zj, 0 < j < n. Podstawiając 
1 


(z — £o) (£j — 21) -+ (£j — £j—ą)(Tj — Lj+1) `: (Lj — In) 


dj = 
oraz a(x) = (£ — £o) -++ (£ — £n) mamy 


węa) = ala) Y EL, (6.4) 


j=0 T — Tj 
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albo A 
2-0 qj (x) f(z;) 
2-0 q;(7) | 
gdzie q;(1) = a;/(x — 1;). W ostatniej równości wykorzystaliśmy fakt, 


że a(x) = (007-ę4;(%)) 77, co wynika z (6.4) po podstawieniu f = 1. 


w(x) = 


BAZA POTĘGOWA (NATURALNA) 


Znacznie prościej można obliczyć wartość wielomianu interpolacyj- 
nego, (a także jego pochodnych), gdy jest on dany w najczęściej uży- 
wanej bazie potęgowej, p;(x) = x”, Vj. Jeśli bowiem 


wple) = ao + at +: + ana”, 


to również 


wple) = (1: (ant + an-1)£ + an-2)£ +: + a1)£ + a0, 


co sugeruje zastosowanie następującego schematu Hornera do obliczenia 
w(x): 


Un := An; 
for j := n—1 downto 0 do 
Uj [= Uj+1 *T + Gy, 


Po wykonaniu tego algorytmu wy(x) = vo. Schemat Hornera wymaga 
wykonania tylko n mnożeń i n dodawań. Ma on również głębszy sens, 
zob. Ćw. 6.2. 

Zauważmy jednak, że w przypadku bazy potęgowej macierz (2; );,—0 
układu (6.1) jest pełna. Jest to tzw. macierz Vandermonde a. Obli- 
czenie współczynników wielomianu interpolacyjnego w bazie potęgowej 
bezpośrednio z tego układu, stosując jedną ze znanych nam już metod, 
kosztowałoby rzędu n operacji arytmetycznych. 


BAZA NEWTONA 

Rozwiązaniem pośrednim, które łączy prostotę obliczenia współ- 
czynników z prostotą obliczenia wartości w(x) jest wybór bazy New- 
tona, 
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W tym przypadku współczynniki rozwinięcia wielomianu interpolacyj- 
nego będziemy oznaczać przez bj, 


wf = M. bjp. 


j=0 


Zwróćmy od razu uwagę na ważną własność bazy Newtona. Jeśli wę, € 
II, jest wielomianem interpolacyjnym dla funkcji f opartym na węzłach 
Zo, X1,...,1j, 0 ŚJ <n, to wfo = bo oraz 


Wjj = Węj-1 + bjPj, 1<j<n. (6.5) 


Wartość w;(z) można obliczyć stosując prostą modyfikację algo- 
rytmu Hornera: 


WED: 
for j := n—1 downto 0 do 


vj = Ujqą * (x — Ty) + Dy. 


Ponadto układ równań (6.1) jest trójkątny dolny (tzn. p,(x,) = 0 dla 


i < j), co pozwala na obliczenie b, algorytmem z Rozdziału 3.1 o koszcie 


n?, 


f (x0); 


Ts a =l F n do 
Se= żz0bspslzy)) e): 


Zwykle jednak do obliczania współczynników b; stosuje się nieco 
inny algorytm, który teraz predstawimy. 


6.3 Algorytm różnic dzielonych 


Różnicę dzieloną funkcji f opartą na różnych węzłach to, ty, ... „ts, gdzie 
s > 1, definiuje się indukcyjnie jako 


f(ti,t2,...,ts) e f(to,ti,...,ts—1) 
t, — to i 


f(to,tu,...,ts) = (6.6) 


Zachodzi następujące ważne twierdzenie. 
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Twierdzenie 6.1 Współczynniki b; wielomianu interpolacyjnego New- 
tona dla danej funkcji f dane są przez różnice dzielone f w węzłach 
Lo, Li,- .-, Lj, t2N. 


b; = f(£o, £1,- .., £3), OSJ <n. 


Dowód Dla 0 <i<j<n, oznaczmy przez w,, wielomian z I;—; in- 
terpolujący f w węzłach zi, £i+1,.-.-, £j. Wtedy ma miejsce następująca 
równość (i < j): 
— TOW; — (z — z;Jwi; l£ 
w; (z) = (x xi)w +15(7) (x zj) J :( ) Yr. (6.7) 


Tj — Ti 


Aby ją pokazać wystarczy, że prawa strona tej równości, którą ozna- 
czymy przez v(x), przyjmuje wartości f(x.) dla £ = £„ i < s <j. 
Rzeczywiście, jeśli ¿+ 1 < s < j — 1 to 


(1, — Ti) f (£5) — (£s — T4) f (£s) 


Tj — Ti 


v(z,) = = f(zs). 


Ponadto 
ESC SEE poz) 


dj — ti 
oraz podobnie v(x;) = f(x;). Stąd v jest wielominem z II,_; interpolu- 
jącym f w węzłach z,, i < s < j, czyli Wij = v. 
Dalej postępujemy indukcyjnie ze względu na stopień n wielomianu 
interpolacyjnego. Dla n = 0 mamy oczywiście bo = f (£o). Niech n > 1. 
Ponieważ, jak łatwo zauważyć, 


WOn(T) = Wom-1(£) + bnpn(z), 


z założenia indukcyjnego mamy b; = f(z£o,..., £j) da 0< j <n-l. 
Aby pokazać podobną równość dla b„, wykorzystamy równość (6.7) z 
i =0ij =n, która ma wtedy postać 


(£ — zo) W nlx) — (z — En)Won-i (2) 


Won(T) = 
Ln — To 
Zauważmy teraz, że b, jest współczynnikiem przy x” w wielomianie 


won. Z założenia indukcyjnego wynika, że współczynniki przy z”! w 
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wielomianach Win i Won-1 Są ilorazami różnicowymi opartymi odpo- 


wiednio na węzłach 11,...,Tqn 1 £0,...,Tn_1. Stąd 
f (£1, ---,En) — f (£0, ---,En-1) 
bn = - — - — = I dode): 
dn WO 


co kończy dowód. 


Różnicę dzieloną f(£o, T1,...,%,) można łatwo obliczyć na podsta- 
wie wartości f(x;), 0 < j < n, budując następującą tabelkę: 
zo f(z0) 


xı f(z1) f(%0,11) 
1a f(%a) f(£1,£2)  f(%o,T1, £2) 


En f(n) FUn=nZn) f(En-2,En-1;, En) > J(Z071,:-:,Tn). 
Zauważmy przy tym, że “po drodze” obliczamy f(£;, £i+1,.-., £j) dla 
wszystkich 0 < i < j < n, a więc w szczególności również interesu- 
jące nas różnice dzielone f(£o,£%1,..., xj). Stąd i z Twierdzenia 6.1 na- 
tychmiast wynika algorytm obliczania współczynników b, wielomianu 
interpolacyjnego Newtona. Po wykonaniu następującego programu, 

for j:=0 to n do bjj] = f(x[j]); 
for j := 1 to n do 
for k := n downto j do 


bik] := (blk] — blk — 1))/(zlk] — z|k — j]), 


współczynniki b, zostaną umieszczone w tablicy bļ[j]. 


6.4 Przypadek węzłów wielokrotnych 


Uogólnieniem rozpatrzonego zadania interpolacji jest zadanie interpo- 
lacji Hermite'a. Zakładamy, że oprócz (różnych) węzłów x, dane są 
również ich krotności n,, 0 < j < k, przy czym 5% nj =n+1. Należy 
skonstruować wielomian wy € II, taki, że 


WPa j= O) da O0<i<nm-10<j<k. 


Oczywiście zakładamy przy tym, że odpowiednie pochodne funkcji f 
istnieją. 
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Lemat 6.2 Zadanie interpolacji Hermite’a ma jednoznaczne rozwiąza- 
nie. 


Dowód Istnienie i jednoznaczność rozwiązania można uzasadnić tak 
samo jak w przypadku węzłów jednokrotnych. Przedstawiając wielo- 
mian w dowolnej bazie otrzymujemy układ n + 1 równań z n+ 1 nie- 
wiadomymi, który dla zerowej prawej strony ma jedynie rozwiązanie 
zerowe. Inaczej bowiem istniałby wielomian niezerowy stopnia nie więk- 
szego niż n, który miałby zera o łącznej krotności większej niż n. © 


Nas oczywiście interesuje konstrukcja wielomianu wę. W tym celu 
ustawimy węzły z, w ciąg 


(UGUls sam ICE (Dose T io fu Lo) 
——— ——>;-— 
no nı Nk 

i zdefiniujemy uogólnioną bazę Newtona w II, jako 

polz) = 1, 

pj(z) = (x -— To)(x— T1): (£ — g;j-1), ISJ <n. 

Uogólnimy również pojęcie różnicy dzielonej na węzły powtarzające 

się kładąc 


A+(2 Za fO- (z;) 
J (To bia dy) = Gz 
Ti ROC = TERRE B= 
POSTA REGAT A= f (Ziyi j) F j-1) 
Tj — Ti 
dla z; Æ t;. Zauważmy, że przy tej definicji różnice f(Z;,...,£;) mo- 


żemy łatwo obliczyć stosując schemat podobny do tego z przypadku 
węzłów jednokrotnych. 


Twierdzenie 6.2 Współczynniki b; wielomianu interpolacyjnego Her- 
mite’a w bazie Newtona, 


wj) = X bpi); 


dane są przez odpowiednie różnice dzielone, tzn. 


bkRJlvzessT e, DJ 
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Dowód Dowód przeprowadzimy podobnie jak dla węzłów jednokrot- 
nych. Niech w,, € II,_, oznacza wielomian interpolacyjny Hermite'a 
oparty na (być może powtarzających się) węzłach Zi, T;44,...,Tj. To 
znaczy, wW; j interpoluje f w węzłach x, takich, że z, występuje w ciągu 
Ti,- .. Tj, a jego krotność jest liczbą powtórzeń z, w tym ciągu. 
Zauważmy najpierw, że dla z; Æ £; zachodzi znany nam już wzór 
(6.7), 
Wo TE kam IE (6.8) 


Tj — Ti 


Rzeczywiście, oznaczmy przez v(x) prawą stronę powyższej równości. 
Dla k mniejszego od krotności danego węzła z, w ciągu 7;,... Zj, mamy 


AS (żzj= i (zs), a ponieważ 


k (wir? (£) — wizy (0) 


(k) — 
BBE Ż;— Dy 
= k = k 
. 6- zw) (a) — (z — żywy) (2) 
| Ż;— 
i = (k) z (k) 
Wuja (Es — Ti)Wwiti (Ts) — (Vs — Tj)wij-1(2s) 
īj— 


Korzystając z tego wzoru sprawdzamy, że v spełnia odpowiednie wa- 
runki interpolacyjne, a stąd w; j = v. 

Dalej postępujemy indukcyjnie ze względu na n. Dla n = 0 mamy 
bo = f (£o). Dla n > 1 wystarczy pokazać, że b, = f(Zo,71,...,£,). W 
tym celu rozpatrzymy dwa przypadki. 

Jeśli to = Zn to mamy jeden węzeł zo o krotności n + 1. Wielomian 
interpolacyjny jest wtedy postaci 


n f(r | 
wę(z) => t o) (x — xto), 
j=0 F 
a stąd b, = f©(zg)/(n!) = f(£o,..., £o). Jeśli zaś Zo £ Z; to równość 
1 
n+ 
bn = f(zo,T1,..., Zn) wynika z (6.8) po podstawieniu i = 0i j =n, 


oraz z założenia indukcyjnego. 


74 ROZDZIAŁ 6. INTERPOLACJA WIELOMIANOWA 


Uwagi i uzupełnienia 


U. 6.1 Algorytm Hornera okazuje się optymalny. Każdy inny algorytm ob- 
liczający dokładną wartość wielomianu znając jego współczynniki wymaga 
wykonania co najmniej n mnożeń i n dodawań. Algorytm Hornera jest też 
numerycznie poprawny, zob. Ćw. 6.1. 


U. 6.2 Mając obliczone współczynniki b; wielomianu interpolacyjnego w 
postaci Newtona, można łatwo przejść do jego współczynników a; w bazie 
potęgowej. Algorytm wynika bowiem z następującego wzoru rekurencyjnego. 

Dlai=n,n—1,...,0, niech > i < j < n, będą odpowiednimi współ- 
czynnikami w rozwinięciu 


b; +b; 1 (©— zi) +- -+brn(x— zi) (£—£n-1) = a +a() ES . kati, 


Wtedy ać = b, oraz 


+1 i+1 EARE 
al Sa | = graj, iLj<n-lL. 
Uwzględniając fakt, że poszukiwanymi współczynnikami są a; = a. 0< 
j < n, algorytm może być zrealizowany następująco: 


for j:=0 to n do ajj] := b[j]; 
for i := n—1 downto 0 do 
for j :=i to n—1 do 

a|j] := alj] — zļ[i] *a|j +1]. 


U. 6.3 Okazuje się, że przy realizacji w fi, algorytmu różnic dzielonych 
istotną rolę odgrywa porządek węzłów. Można pokazać, że algorytm liczenia 
f(to,...,tn) jest numerycznie poprawny ze względu na dane interpolacyjne 
f O(t), o ile węzły są uporządkowane nierosnąco lub niemalejąco, zob. Ćw. 
6.5. 


U. 6.4 Zauważmy, ze pojęcie różnicy dzielonej formalnie zdefiniowaliśmy je- 
dynie dla ciągu węzłów postaci 10,...,T70,01,..-,71,...,Tk,..., Th, gdzie Tj 
są parami różne. Tą definicję można rozszerzyć do dowolnego ciągu węzłów. 
Można bowiem powiedzieć, że f(to,t1,...,tn) jest współczynnikiem przy z” 
wielomianu wy,...«„ € In interpolującego f w węzłach t; (uwzględniając 
krotności). Równoważnie, 

wy) 


t0--utn 
ROEE eb) =“ i . 
n: 
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Cwiczenia 
Ćw. 6.1 Pokazać numeryczną poprawność algorytmu Hornera obliczania 


wartości wielomianu ze względu na dane współczynniki a; tego wielomianu. 


Ćw. 6.2 Pokazać, że algorytm Hornera obliczania wartości w(ć) wielomianu 
danego w postaci potęgowej jest jednocześnie algorytmem dzielenia tego wie- 
lomianu przez jednomian (x — €). Dokładniej, jeśli w(x) = © ?_o aga) to 


(z) = (Yua) +(r-6) + w, 
j= 


gdzie vj są zdefiniowane tak jak w algorytmie Hornera. 


Ćw. 6.3 Zaproponować możliwie tani algorytm obliczający współczynniki 
wielomianu w bazie Newtona, dysponując współczynnikami tego wielomianu 
w bazie potęgowej. 


Ćw. 6.4 Pokazać, że dla różnych węzłów tg,...,tn mamy 
n 
F(tj) 
f(to,t,...;tn) = 
e 2 (tj — to) + (tj — tj-1) (tj — tj+1) (tj — tn) 
Ćw. 6.5 Niech ty < ti < ... < tn. Pokazać, że realizując w fi, algo- 
rytm różnic dzielonych, zamiast dokładnej wartości f(to,t1,...,tn) otrzymu- 


jemy fly (Flto, tes o, która jest dokładną różnicą dzieloną dla danych 
f(t;)(1 + ez), gdzie |e;| < Krv, 0< J<n. 

Wskazówka Rozpatrzyć najpierw proste przypadki n = 1, 2. Dla dowolnego 
n skorzystać z Ćw. 6.4. 


Ćw. 6.6 Uzasadnić, że iloraz różnicowy jest funkcją symetryczną, tzn. dla 
dowolnej permutacji tio,- --, tis węzłów łg,...,£s mamy 


UN tir) = J tozenta): 


Ćw. 6.7 Niech w będzie wielomianem stopnia dokładnie n. Pokazać, że dla 
ustalonych ty,to,...,ts funkcja 


W... ts (t) = w(t, ti, t2, ... sta) tE R, 


jest wielomianem stopnia dokładnie n — s dla s < n, oraz jest zerem dla 
s>n+l. 


16 ROZDZIAŁ 6. INTERPOLACJA WIELOMIANOWA 


Ćw. 6.8 Pokazać, że jeśli funkcja f : D + R jest n-krotnie różniczkowalna 
na D w sposób ciągły, f € C™(D), to 


f(to, ti, mięty to) = BR f(to,t1, sasz Uh ZAJE 


Ćw. 6.9 Stosując algorytm różnic dzielonym znależć wielomian w postaci 
Hermite'a interpolujący funkcję f(x) = xt w trzech dwukrotnych węzłach 
równych 0,1i 2. 


CE 


Rozdział 7 


Interpolacja a aproksymacja 
funkcji 


Gdy mamy do czynienia z funkcją, która jest “skomplikowana” to często 
dobrze jest zastąpić ją funkcją “prostszą”. Mówimy wtedy o aproksyma- 
cji (przybliżaniu) funkcji. Funkcję musimy również aproksymać wtedy, 
gdy nie jesteśmy w stanie uzyskać pełnej o niej informacji. Na przykład, 
gdy funkcja reprezentuje pewien proces fizyczny to często zdarza się, 
że dysponujemy jedynie ciągiem próbek, czyli wartościami tej funkcji 
w pewnych punktach. Jasne jest, że chcielibyśmy przy tym, aby błąd 
aproksymacji był możliwie mały. 

Z tego punktu widzenia, intepolacja wielomianowa może być trak- 
towana jako jeden ze sposobów aproksymacji funkcji, opartym na prób- 
kowaniu. Naturalnym staje się więc pytanie o błąd takiej aproksymacji. 


7.1 Błąd interpolacji wielomianowej 


Niech zo, £1,..., £n będą (niekoniecznie różnymi) węzłami należącymi 
do pewnego (być może nieskończonego) przedziału D C R. Dla da- 
nej funkcji f : D — R, przez wę oznaczymy, tak jak w poprzednim 
rozdziale, wielomian interpolacyjny stopnia co najwyżej n interpolu- 
jący f w zadanych węzłach. W przypadku węzłów wielokrotnych jest 
to oczywiście wielomian interpolacyjny Hermite'a. 


TT 
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Lemat 7.1 Dla dowolnego punktu z € D błąd interpolacji w £ wyraża 
się wzorem 


JE) — w;(2) = (Z— z0)(Z—14):*: (Z— 1n)f(%o,T4,...,En,7). (7.1) 
Jeśli ponadto f € C"+V(D), czyli pochodna f+) w D istnieje i jest 
ciągła, to 
Fr) 

(n+ 1)!’ 


gdzie £ = £(1) jest pewnym punktem należącym do najmniejszego prze- 


działu zawierającego punkty £o, £1, ..., Zn, T. 


f(2) — w;(i) = (I — £o) (T — x1): (E — £n) 


Dowód Możemy założyć, że z nie jest żadnym z węzłów z,,0<j< 
n. Niech w; € Iin} będzie wielomianem interpolacyjnym funkcji f 
opartym na węzłach 1tg,...,%, i dodatkowo na węźle z. Mamy wtedy 


Wę(z) = wple) + (x — zo) (x — x1) -++ (£ — 2,)f(X0,£1,.-.;Ty,Z), 


a ponieważ z warunku interpolacyjnego f(z) = wy(z), to mamy też 
(7.1). 
Aby pokazać drugą część lematu, rozpatrzmy funkcję % : D > R, 


v(z) = f(x) — wę(z) 


= f(x) — wę(z) — (x — 20)(z — 14): (£ — En) f (20, - -, 2n, 2). 


Z warunków interpolacyjnych na wy € II,,, wynika, że funkcja V ma 
punkty zerowe o łącznej krotności co najmniej n + 2. Wykorzystując 
twierdzenie Rolle'a wnioskujemy stąd, że y’ ma zera o łącznej krotności 
co najmniej n+ 1, Y” ma zera o łącznej krotności co najmniej n, itd. W 
końcu funkcja ("1 zeruje się w co najmniej jednym punkcie £ = £(z) 
należącym do najmniejszego przedziału zawierającego t0,11,...,dn,T. 
Wobec tego, że u = 0, a (n + 1)-sza pochodna wielomianu (z — 


©o):*:(1— £n) wynosi (n + 1)!, mamy 
0 = YDE) = SOE) — (n + 1)! (T0, , 2n, 2). 
Stąd 


(n+1) 
ftot een amna) = R 
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co kończy dowód. 


Zwykle interesuje nas nie tyle błąd w ustalonym punkcie z € D, ale 
na całym przedziale D. Zakładając teraz, że przedział D jest domknięty, 
czyli 

D = ja,b] 
dla pewnych —oo < a < b < +w, błąd ten będziemy mierzyć w normie 


jednostajnej (Czebyszewa). Dla funkcji ciągłej g : [a,b] + R, norma ta 
jest zdefiniowana jako 


lallen) = maz|g(0)|. 
Niech Fy([a,b]), gdzie r > 0, będzie klasą funkcji 
Fyrlla,t) = {f € Clad) : | lleqawy < M}, 
gdzie 0 < M < oo. Mamy następujące twiedzenie. 


Twierdzenie 7.1 Załóżmy, że każdą funkcję f € EFqx,(|a,b|) aproksy- 
mujemy jej wielomianem interpolacyjnym wę € Il, opartym nar+1 
węzłach 10,...,«, € [a,b]. Wtedy maksymalny błąd takiej aproksymacji 
wynosi 


e(Fy(la,b]); £o, £1,- - , £r) 


max — w à 

PE) IF — wz||C(a,t)) 
M 

"Go 22M EZ 


Dowód Oszacowanie górne wynika bezpośrednio z Lematu 7.1, bo- 
wiem dla f € Fy([a,b|) mamy 


If = wy||cfat) = max |F(z) — wę(z)| 
(r+1) 
= ae (x = zo) SA (x — £r )| If 7 e 


Ge DE e-e) (e-e) 


S0ROÓZDZIAŁ 7. INTERPOLACJA A APROKSYMACJA FUNKCJI 


Z drugiej strony zauważmy, że dla wielomianu v(x) = Ma"*'/(r+1)! 
mamy v € Fz,([a,b|) oraz 


M SBIEGAC "RED a 


co kończy dowód. 


Zauważmy, że błąd aproksymacji e(F5,([a,b)); £o, . .-, £r) w istotny 
sposób zależy od wyboru węzłów xj. Naturalne jest więc teraz następu- 
jące pytanie. W których punktach x; przedziału [a,b] należy obliczać 
wartości funkcji, aby błąd był minimalny? Problem ten sprowadza się 
oczywiście do minimalizacji wielkości maxy<q< |(V — xo): +- (£ — 2,)| 
względem węzłów xj. 


Twierdzenie 7.2 Błąd aproksymacji Fy(|a,b])(zo,***,«,) jest mini- 
malny gdy węzły 


b — 23 +1 b 
m= = cos (5 r) Si j O<j<r. 
2r+2 


Ponadto, dla węzłów optymalnych x; mamy 


Fpa a = Zz (5) 


Dowód tego twierdzenia opiera się na własnościach pewnego waż- 
nego ciągu wielomianów, który teraz przedstawimy. 


7.2 Wielomiany Czebyszewa 


Ciąg {Tk}k>0 wielomianów Czebyszewa (pierwszego rodzaju) zdefinio- 
wany jest indukcyjnie jako 


T(z) = 1, 
Tı(£) = m 


Tra) dla k>1. 


I 
N 
8 
= 

> 

l 
= 

L 

R 
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Zauważmy, że Th jest wielomianem stopnia dokładnie k o współ- 
czynniku przy x” równym 2*7! (k > 1). Ponadto wielomian Ty można 
dla |x| < 1 przedstawić w postaci 


T(x) = cos(k arccos z). (7.2) 


Rzeczywiście, łatwo sprawdzić, że jest to prawdą dla k = 0,1. Stosując 
podstawienie cost = z, 0 < t < m, oraz wzór na sumę cosinusów 
otrzymujemy dla k > 1 


cos((k + 1)t) = 2- costcos(kt) — cos((k — 1)t), 


co jest równoważne formule rekurencyjnej dla 754,1. 
Ze wzoru (7.2) wynikają również inne ważne własności wielomianów 
Czebyszewa. Norma wielomianu Czebyszewa na |—1, 1] wynosi 


|Fzlcq—11y 5 max |74(x)| = 1 


—1<x<1 


i jest osiągana w k + 1 punktach tego przedziału równych 


Y = cos (r), G<jE, (7.3) 


przy czym Ty(y;) = (—1)/. 
W końcu, k-ty wielomian Czebyszewa 14 ma dokładnie k pojedyn- 
czych zer w [—1, 1] równych 


23 +1 
z = cos (r), 0<J<k-1. 


Konsekwencją wymienionych własności jest następujący lemat, któ- 
ry jest kluczem do dowodu Twierdzenia 7.2. 

Przez II, oznaczymy klasę wielomianów stopnia k o współczynniku 
wiodącym równym 1, tzn. 


M, = {w E Mp: wr) =+}. 


Lemat 7.2 Niech k > 1. W klasie Il, minimalną normę jednostajną 
na przedziale [—1,1] ma wielomian w* = 27 *TĘ, tzn. 


. + 1 
P llwlleq-1,1) = ||w |C[-11) T 9k-1' 
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Dowód Zauważmy najpierw, że w* € II; oraz |w*|ef—11p = 27%. 
Wystarczy więc pokazać, że norma każdego wielomianu z II, jest nie 
mniejsza niż 217%, 

Załóżmy, że jest przeciwnie, tzn. istnieje wielomian w € II, taki, że 


1 * 
lwlleq-1p < z w lleq-1,1- 
Rozpatrzmy funkcję Y = w* — w. Ponieważ dla punktów “maksy- 
malnych” zdefiniowanych w (7.3) mamy w*(y-;) = (—1)72'7* oraz 


[w(yr-;)| < 277%, to 


>0 j-parzyste, 
<0 J-nieparzyste. 


U(Ykj) 


0< j <k. Stąd v ma co najmniej jedno zero w każdym z przedziałów 
(Yi, Yi+1) dla 0 < i < k — 1, czyli w sumie k zer. Z drugiej strony, % jest 
wielomianem stopnia co najwyżej k — 1 (bo współczynniki przy z* w 
wielomianach w* i w redukują się), a więc Y = 0 i w* = w. 


Dowód Twierdzenia 7.2 wynika teraz bezpośrednio z Lematu 7.2. Za- 
uważmy bowiem, że wielomian (x — tg)(z — 14): -- (1 — zr) jest w klasie 
I+. Stąd dla [a,b| = [1,1] optymalnymi węzłami są zera z; wielo- 
mianu Czebyszewa, przy których 


dą (x) f 


(x — zo) (x — z1) (£ — 2,) = J 


Jeśli przedział [a,b] jest inny niż |—1, 1], należy dokonać liniowej za- 
miany zmiennych tak, aby przeszedł on na [—1, 1]. Bezpośrednie spraw- 
dzenie pokazuje, że w klasie II.,, minimalną normę Czebyszewa na 
przedziale [a,b] ma wielomian 


b—asr+1 „(2x — (a+b) 
z) w (ZJ. 


(2) = ( 


Stąd 


Z 1 — Te 


[wa oll Cia) = > ży 


i węzły x} są optymalne. © 
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7.3  Interpolacja kawałkami wielomiano- 
wa 


W praktyce interesuje nas konstrukcja nie jednej, a raczej ciągu aprok- 
symacji, które przybliżałyby daną funkcję f z coraz większą dokładno- 
ścią. Cel ten można spróbować osiągnąc przez wzięcie jako n-tej aprok- 
symacji wielomianu interpolacyjnego wę, € II, opartego na pewnych 
węzłach z;, 0 < j < n. Rzeczywiście, jeśli f jest na [a,b] nieskończenie 
wiele razy różniczkowalna to n-ty błąd dla dowolnego układu węzłów 
jest ograniczony przez 


|FET||cdat) 


GE DI (b — a)", (7.4) 


Ff — Wgn||C(aż) < 


a więc szybko dąży zera wraz z n + œ, gdy np. wszystkie pochodne 
funkcji f są wspólnie ograniczone przez pewną stałą. Tak jest np. dla 
f(x) = e”. Jeszcze lepsze oszacowanie dostaniemy, gdy wybierzemy 
węzły Czebyszewa. 

Podejście to ma jednak pewne wady. Zwiększając stopień wielo- 
mianu interpolacyjnego sprawiamy, że aproksymacja, a także obliczenie 
współczynników wielomianu interpolacyjnego, stają się coraz bardziej 
skomplikowane, a tego chcielibyśmy uniknąć. Poza tym, założenie o 
istnieniu i wspólnej ograniczoności wszystkich pochodnych funkcji nie 
zawsze jest spełnione (np. w rozpatrywanej już klasie F5,(|a,b])), co po- 
woduje, że oszacowanie (7.4) nie zawsze jest prawdziwe, a nawet błąd 
If = wznllcav) Może nie maleć do zera, zob. U. 7.1. 


Podamy teraz inny sposób konstrukcji ciągu aproksymacji, który 
zapewnia zbieżność w klasie Fx,(|a,b|), ale nie wymaga nadmiernego 
zwiększania stopnia wielomianu interpolacyjnego. 


Postąpimy zgodnie ze znaną zasadą “dziel i rządź”. Podzielimy 
przedział [a,b] na k równych podprzedziałów [t;_1,t,], 1 < i < k, i 
zastosujemy na każdym z nich interpolację wielomianową stopnia r. 
Dokładniej, dla każdego i wybieramy (być może wielokrotne) węzły 
Tij E€ [yt], 0 < j < r, po czym aproksymujemy f na [t, 4, £| jej 
wielomianem interpolacyjnym opartym na tych węzłach. Tak skonstru- 
owaną kawałkami wielomianową aproksymację funkcji f : [a,b] > R 
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oznaczymy przez wę. Zauważmy, że wę wykorzystuje na całym [a,b] 
węzły o łącznej krotności n spełniającej nierówność 


n śk(r+1). (7.5) 


Twierdzenie 7.3 Błąd aproksymacji kawatkami wielomianowej w kla- 
sie funkcji Fxy,(|a,b|) jest ograniczony przez 


M b—asr+1 lr+1 
=le gs LEANE 
JeF? llab) If Wyę||C(l d]) = FED k ) = > 


gdzie 
„7, M (b = a)" ti(r + r 
(r +1)! 
Jeśli ponadto węzły xij, 0 < j <r, są na każdym przedziale |t;_1, t:l 
węzłami Czebyszewa, to dla odpowiedniej aproksymacji wę mamy 


max IF — wple = e 3 ZA Ey” 


feFr, (la (r+ 1)! 4k n 
gdzie 
g Si _ 2M(b— at (r +1)" 
A2 4"+1(r + 1)! 


Dowód Niech f € Fy(|a,b]) i z € [ti-1, t:i]. Z Lematu 7.1 wynika, że 


EERE] 


G+ DI (£ — zio) (£ — £11) (£ — Eir), 


|f) — wę(z)| = 


a ponieważ |x — 1; ;| < ti — ty 4 = (b — a)/k, to 


PORTOT N 


(r+1)!\ k 
Ta nierówność wraz z (7.5) dowodzi pierwszą część twierdzenia. 


Dla dowodu drugiej części wystarczy wykorzystać Twierdzenie 7.2, 
aby zauważyć, że 


IMCZEKELORCJE 
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oraz, że dla f(x) = Mz'*"/(r + 1)! mamy powyżej równość. 


Otrzymaliśmy, że interplacja kawałkami wielomianami stopnia r ko- 
rzystająca z n wartości funkcji daje w klasie funkcji Fxy,(|a,b|) błąd 
rzędu n""*V. Błąd ten zbiega więc do zera, gdy liczba n wartości funk- 
cji rośnie do nieskończoności, przy czym charakter zbieżności (tzn. przy 
zaniedbaniu stałej przy n""*V) nie zależy ani od długości przedziału 
la, b], ani od doboru *wewnętrznyh” węzłów zx,,,0<j<r. 

Można pokazać, że powyższa aproksymacja jest w klasie Fxy,([a,b]) 
optymalna, tzn. błąd dowolnej innej metody aproksymacji korzystającej 
z wartości funkcji w n punktach dąży do zera nie szybciej niż n""*V; 
zob. U.7.2. 


Uwagi i uzupełnienia 


U. 7.1 Niech f(x) = (1+a27)"1. Okazuje się, że jeśli wy, jest wielomia- 
nem interpolacyjnym dla f opartym na węzłach równoodległych przedziału 
[-5,5], tzn. z; = —5 + 10j/n, 0 < j < n, to ciąg wę„(z) jest zbieżny do 
f(x) dla |x| < 3.68... i rozbieżny dla |x| > 3.68.... To znaczy, błąd interpo- 
lacji funkcji f na całym przedziale [—5,5) nie maleje do zera, gdy stopień 
wielomianu interpolacyjnego rośnie do oo. 


U. 7.2 Pokażemy teraz, że aproksymacja funkcjami kawałkami wielomia- 
nowymi jest optymalna w F;y([a,b|) ze względu na rząd zbieżności błędu. 
Dokładniej pokażemy, że dla dowolnej aproksymacji A(f) funkcji f wykorzy- 
stującej wartości f lub jej pochodnych w co najwyżej n różnych punktach, 
tzn. dla A postaci 


A(f) = (Ff(z1), 0% fo (m), ea (En), 253 fTFV(2,)) 


mamy 
sup |f- A(Ollogary Z CU MnC, (7.6) 
feFy, (la,b]) 
gdzie C1 > 0 jest pewną stałą niezależną od A, M in. 
W tym celu wybierzmy dowolną niezerową funkcję v : R — R spełnia- 
jącą warunki: 


1. y € Fi (R), 
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2 y(x) = 0 dla z ¢ [0,1]. 


Niech s będzie takim indeksem, że hs = £s4+1 — z, > (b— a)/(n +1). Wtedy 
dla funkcji Y;, i = 1,2, zdefiniowanych jako 


(e) = vla) = M(E) 


sS 
mamy p; € Fy([a,b]) oraz informacja o y; jest zerowa, yP (z) = 0 dla 
0<p<r+1,1<j< n. Stąd go = 9(0,...,0) jest aproksymacją dla obu 
<< — 
n(r+2) 
funkcji y; i z nierówności trójkąta dla norm dostajemy 


max ||%1 — $o||E(tat)> 192 — poll Cab} 
> lyi — %ellC([at/2 = IllC(at) 


b—ayr+1 
r+1 
= Mh" ly|c(o) Z ME ep Co.) 


Nierówność (7.6) zachodzi więc z C14 = (2(b — aj)" +"|ellc(0.1)- (Zob. też 
Ćw. 8.3.) 


U. 7.3 Załóżmy, że funkcje z klasy Fy,([a,b|) chcemy aproksymować z do- 
kładnością e > 0 w normie Czebyszewa używając możliwie mało obliczeń 
wartości funkcji. Z Twierdzenia 7.3 i U. 7.2 wynika, że musimy wtedy obli- 
czać co najmniej rzędu (1/e)7*! wartości funkcji. Co więcej, aproksymacja 
kawałkami wielomianowa daje błąd e obliczając właśnie tyle wartości funkcji. 
Mówimy, że złożoność informacyjna zadania jest rzędu (1/e)"*". 


Cwiczenia 


Ćw. 7.1 Niech dane będą € > 0 i funkcja f : [a,b] + R, która jest nieskoń- 
czenie wiele razy różniczkowalna i wszystkie jej pochodne są na [a,b] ograni- 
czone przez M. Napisać program znajdujący stopień n oraz współczynniki 
w bazie Newtona wielomianu wę, € II, interpolacyjnego f z błędem 


IF = ©gn|C(ab]) < E- 


Rozpatrzyć dwa przypadki: gdy węzły interpolacji są równoodległe, oraz gdy 
węzły są czebyszewowskie. 
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Ćw. 7.2 Pokazać, że dla k > I mamy 


ENS E 


2 2 


Ty(z) = 
gdzie TĄ jest k-tym wielomianem Czebyszewa. 


Ćw. 7.3 Pokazać, że wielomiany Czebyszewa tworzą układ ortogonalny w 
przestrzeni LQ(—1,1), gdzie waga p(x) = (1 — 27)71/2, tzn. iloczyn skalarny 


yo (BABE E a r 
(Ti, Tj) = Je der dx = i o 0, 


Ćw. 7.4 Wielomiany Czebyszewa 17%, 0 < k < n, tworzą bazę w przestrzeni 
II. Podać algorytm obliczający współczynniki danego wielomianu w bazie 
Czebyszewa na podstawie jesgo współczynników w bazie potęgowej, albo w 
bazie Newtona. 


Ćw. 7.5 Podać algorytm typu algorytmu Hornera oblicząjący wartość wie- 
lomianu w(x) danego przez jego współczynniki w bazie Czebyszewa {Tę |k>0. 


Ćw. 7.6 Niech w będzie wielomianem, którego rozwinięciem w bazie Cze- 
byszewa jest w(x) = X g-o akTk(x). Pokazać, że wielomian 


nl 
w(x) = X axTk(z) 
k=0 


najlepiej przybliża w w normie Czebyszewa na przedziale [—1,1] spośród 
wszystkich wielomianów stopnia co najwyżej n — 1, oraz 


[w w|cq-11) = |danl. 
Czy w(x) = Da aęTy(x) najlepiej przybliża w w tej samej normie w 
klasie Il„_2? 
Ćw. 7.7 Załóżmy, że wielomian w(x) = X p=0 azTk(x) aproksymujemy wie- 
lomianem w(x) = 24-65 a@kTk(x). Uzasadnić, że 


n 


lw- wlleqip < 2, lal. 
k=n-i+1 


Kiedy w tej nierówności mamy równość? 
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Ćw. 7.8 Niech dany będzie przedział skończony [a, b) nie zawierający zera, 
oraz liczba c. Niech 


Ü, = {w € Mp: w(0)=c]. 
Pokazać, że w klasie IL, najmniejszą normę Czebyszewa na [a,b| ma wielo- 
mian kaei 
x—(a+ 
= 
a a 
T, (33) 

i wynosi ona 


lwli) = 


Jakie jest rozwiązanie, gdy 0 € Ja, b]? 


Rozdział 8 


Interpolacja funkcjami 
sklejanymi 


Interpolacja kawałkami wielomianami interpolacyjnymi daje mały błąd 
aproksymacji, ma jednak również pewne wady. Zauważmy, że w punk- 
tach t; podziału przedziału [a,b] na podprzedziały funkcja aproksymu- 
jąca wy nie jest w ogólności ciągła. W najlepszym wypadku, nawet gdy 
węzły x; j na przedziale [t;—1, t:] wybierzemy tak, że £;o = ti-1 Í Zir = ti, 
Vi, to wę(w) co prawda będzie ciągła, ale nie będzie różniczkowalna w 
xi. Ta własność jest niepożądana w przypadku gdy wiemy, że funk- 
cja którą aproksymujemy jest “gładka”. Wtedy lepiej jest zastosować 
innego rodzaju interpolację, np. posługując się funkcjami sklejanymi. 


8.1 Co to są funkcje sklejane? 


W ogólności przez funkcję sklejaną rozumie się każdą funkcję przedzia- 
łami wielomianową. Nas będą jednak interesować szczególne funkcje 
tego typu i dlatego termin funkcje sklejane zarezerwujemy dla funk- 
cji przedziałami wielomianowych i posiadających dodatkowe własności, 
które teraz określimy. 


Niech dany będzie przedział skończony [a,b] i węzły 
a = to < Tı <: < Tn = b, 


przy czym n > 1. 
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Definicja 8.1 Funkcję s : R — R nazywamy funkcją sklejaną rzędu 
r (r > 1) odpowiadającą węzłom xj, 0 < j < n, jeśli spełnione są 
następujące dwa warunki: 


(i) s jest wielomianem stopnia co najwyżej 2r — 1 na każdym z prze- 
działów |xj-1, £5], tzn. Sg, 14] E€ Ilgr=1,1 ŚJ <n, 


(ii) s jest (2r — 2)-krotnie różniczkowalna w sposób ciągły na całej 
prostej, tzn. s € C2 (R). 


Jeśli ponadto 


(iii) s jest wielomianem stopnia co najwyżej r — 1 poza (a,b), tzn. 
S((-00.a]> 5|[b,+oo) € lr—1, 


to s jest naturalną funkcją sklejaną. 


Klasę naturalnych funkcji sklejanych rzędu r opartych na węzłach z; 
będziemy oznaczać przez S,(£o,..., £n), albo po prostu S,, jeśli węzły 
są ustalone. 


Na przykład, funkcją sklejaną rzędu pierwszego (r = 1) jest funkcja 
ciągła i liniowa na poszczególnych przedziałach |1;_1, 1,]. Jest ona na- 
turalna, gdy poza (a,b) jest funkcją stała. Tego typu funkcje nazywamy 
liniowymi funkcjami sklejanymi. 

Najważniejszymi a praktycznego punktu widzenia są jednak funkcje 
sklejane rzędu drugiego odpowiadające r = 2. Są to funkcje, które są 
na R dwa razy różniczkowalne w sposób ciągły, a na każdym z pod- 
przedziałów są wielomianami stopnia co najwyżej trzeciego. W tym 
przypadku mówimy o kubicznych funkcjach sklejanych. Funkcja skle- 
jana kubiczna jest naturalna, gdy poza (a,b) jest wielomianem linio- 
wym. 


8.2 Interpolacja i gładkość 


Pokażemy najpierw ważny lemat, który okaże się kluczem do dowodu 
dalszych twierdzeń. 
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Niech W”(a,b) będzie klasą funkcji f : [a,b] > R takich, że f jest 
(r—1) razy różniczkowalna na [a,b] w sposób ciągły oraz f™ (x) istnieje 
prawie wszędzie na [a,b] i jest całkowalna z kwadratem, tzn. 


W'(a,b) = ff eC"N(fa,b)) : f(x) istnieje p.w. na [a,b] 
oraz f™ € La(a,b) |. 


Oczywiście każda funkcja sklejana rzędu r (niekoniecznie naturalna) 
należy do W” (a,b). 


Lemat 8.1 Niech f € W'(a,b) będzie funkcją zerującą się w węzłach, 
tzn. 
ORA O<S<j<n. 


Wtedy dla dowolnej naturalnej funkcji sklejanej s € Sp mamy 


f fO(x)s"(zx) dz = 0. 


Dowód Dla r = 1 funkcja s' jest przedziałami stała. Oznaczając przez 
a; jej wartość na [x;j—1, xj] dostajemy 


Fastjd = ;f paadz 
I J 


E OSW 


ponieważ f zeruje się w tj. 
Rozpatrzmy teraz Ra r>2. Ponieważ 


(PETE 1) 5 y = fOs (r) a PORE Ar 


r fO(a) z) dz = fe D(z -f pe- )s"+D(z) dz. 


Wobec tego, że s jest poza przedziałem (a,b) wielomianem stopnia co 
najwyżej r — 1 oraz s”) jest ciągła na R, mamy s™ (a) = 0 = s") (b), a 
stąd 


= 0. 


a 


J(e) (w) 
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Postępując podobnie, tzn. całkując przez części r — 1 razy otrzymujemy 
w końcu 


f Worde e fi "Pee de 


Funkcja s") jest przedziałami stała, a więc możemy teraz zastosować 
ten sam argument jak dla r = 1, aby pokazać, że ostatnia całka jest 


równa zeru. 


Funkcje sklejane chcielibyśmy zastosować do interpolacji funkcji. 
Ważne jest więc, aby odpowiednie zadanie interpolacyjne miało jed- 
noznaczne rozwiązanie. 


Twierdzenie 8.1 Jeślin+1 > r, to dla dowolnej funkcji f : [a,b] > R 
istnieje dokładnie jedna naturalna funkcja sklejana s; E€ S, interpolu- 
jąca f w węzłach xj, tzn. taka, że 


sę(z;) = f(z), (ejźm. 


Dowód Pokażemy najpierw, że jedyną naturalną funkcją sklejaną in- 
terpolującą dane zerowe jest funkcja zerowa. Rzeczywiście, jeśli s(1;) = 
Odla0<j<n, to podstawiając w Lemacie 8.1 f = s otrzymujemy 


[ CAO dx = Q. 


Stąd sl”) jest funkcją zerową, a więc s jest wielomianem stopnia co 
najwyżej r — 1 zerującym się w co najmniej n + 1 punktach z,. Wobec 
tego, że n + 1 > r— 1, otrzymujemy s = 0. 

Zauważmy teraz, że problem znalezienia naturalnej funkcji skleja- 
nej s interpolującej f można sprowadzić do rozwiązania kwadratowego 
układu równań liniowych. Na każdym przedziale (1, ,z;], 1 <í <n, 


jest ona postaci 
2r—1 


s(x) = wi(z) = > aijt, 


dla pewnych współczynników a;; € R, a na (—œ0,a] i [b,oo) mamy 
odpowiednio 


r—1 
s(x) = wols) = D aya 
j=0 
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s(x) = Wn (£) = N anyit. 


Aby wyznaczyć s, musimy więc znaleźć ogółem 2r(n + 1) współczynni- 
ków aij, przy czym są one związane (2r — 1)(n + 1) warunkami jedno- 
rodnymi wynikającymi z gładkości, 

Pæ) = whale) 

dla 0<i<nil < k< 2r-—2, oraz n+ 1 niejednorodnymi warunkami 
interpolacyjnymi, 


w,(z;) = f(z) 
dla 0 < i < n. Otrzymujemy więc układ 2r(n + 1) równań liniowych ze 
względu na 2r(n + 1) niewiadomych a; j. 

Naturalna funkcja sklejana interpolująca f jest wyznaczona jedno- 
znacznie wtedy i tylko wtedy, gdy układ ten ma jednoznaczne roz- 
wiązanie. To zaś zachodzi, gdy zero jest jedynym rozwiązaniem układu 
jednorodnego. Rzeczywiście, układ jednorodny odpowiada zerowym wa- 
runkom interpolacyjnym, przy których, jak pokazaliśmy wcześniej, ze- 
rowa funkcja sklejana (której odpowiada a,; = 0, Vi, j) jest jedynym 
rozwiązaniem zadania interpolacyjnego. 


Naturalne funkcje sklejane możemy więc używać do interpolacji 
funkcji. Pokażemy teraz inną ich własność, która jest powodem dużego 
praktycznego zainteresowania funkcjami sklejanymi. 


Twierdzenie 8.2 Niech f e W'(a,b) i niech sę E€ S, będzie naturalną 
funkcją sklejaną rzędu r interpolującą f w węzłach z;,0< j <n. 


Wtedy 
[ 69) * dz > [ (f (2) (8.1) 


Dowód Jeśli przedstawimy f w postaci f = sę + (f — sę), to 


[ 690) a = [ GP) a+ f (G f- sp) O 


2 | SPU slad. 
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Funkcja f—sy jest w klasie W” (a, b) i zeruje się w węzłach zj, 0< j <n. 
Z Lematu 8.1 wynika więc, że f? > (x)(f —s;)"(x)da = 0, a stąd (8.1). 
o 


Wartość całki [,(f0)(«))?da może być w ogólności uważana za miarę 
gładkości funkcji. Nierówność (8.1) możemy więc zinterpretować w na- 
stępujący sposób. Naturalna funkcja sklejana jest w klasie W” (a, b) naj- 
gładszą funkcją spełniającą dane warunki interpolacyjne w wybranych 
węzłach zj. 

Konsekwencją Lematu 8.1 są również inne aproksymacyjne własno- 
ści naturalnych funkcji sklejanych, zob. U. 8.2. 

Jak już wspomnieliśmy, najczęściej używanymi są kubiczne funkcje 
sklejane. Dlatego rozpatrzymy je oddzielnie. 


8.3  Kubiczne funkcje sklejane 


Jeśli zdecydowaliśmy się na użycie kubicznych funkcji sklejanych to 
powstaje problem wyznaczenia sę € S2 interpolującej daną funkcję f, 
tzn. takiej, że s;(1;) = f(x.) dla0 < i < n. W tym celu, na każdym 
przedziale [x;, v,,,| przedstawimy sę w postaci jej rozwinięcia w szereg 
Taylora w punkcie 2;, 


—7.)2 — g)’ 
sp(æ) e w;(z) = a; + b;(z — zz) ! aż Z ł a U , 
i podamy algorytm obliczania a;, bi, c; d; da 0 <i <n-—1. 

Warunki brzegowe i warunki ciągłości dla sę dają nam wg(zg) = 
0 = wh_;(Xn) oraz W; (xii) = wi (£i), czyli 


o = 0, 
CAR dihi = Ciy, O<i<n-2, 
Cn—1 + dn-1ħhn-1 = 0, 


gdzie h; = 1;,4 — zi. Stąd, przyjmując dodatkowo c, = 0, otrzymujemy 


Ci+1 — Gi 


di = h; , 


ieissa i (8.2) 
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Z warunków ciągłości dla s, dostajemy z kolei 
b; + cihi + dzhż/2 = biz, 0<i<n—-2, 
i uwzględniając (8.2) 
bia = bi + hilii + ci)/2, O0<i<n-2. (8.3) 
Warunki ciągłości sę dają w końcu 
a; + bihi + czhż/2 + dzhż/6 = Qi, 0<i<n-2. (8.4) 


Równania (8.2),(8.3) i (8.4) definiują nam na odcinku [a,b] naturalną 
kubiczną funkcję sklejaną. Ponieważ poszukiwana funkcja sklejana s; 
ma interpolować f, mamy dotatkowych n + 1 warunków interpolacyj- 
nych w;(z;) = f(x;), 0O£i<n—1, oraz w, _1(1,) = f(1,), z których 


ai = f(x), O0<i<n—l1. (8.5) 


Teraz możemy (8.4) przepisać jako 


f(za) = f(t) + bihi + cih? + dzhż /6, 


przy czym wzór ten zachodzi również dla i = n — 1. Po wyrugowaniu 
b, i podstawieniu d; z (8.2), mamy 


b; = J (£i, ©;41) + hi(Ci+1 + 2c;)/6, 0 = 4 < n — 1, (8.6) 


gdzie f(£i, £i+1) jest odpowiednią różnicą dzieloną. Możemy teraz po- 
wyższe wyrażenie na b; podstawić w równaniach (8.3), aby otrzymać 


cihi/6 + ciq1 (hi + hi+1)/3 + Cişihi41/6 EJ) — f (Ti, Ta 
Wprowadzając oznaczenie 
G; = &/6, (8.7) 
możemy to równanie przepisać jako 


hi hi 
ci + 20h + AE — cha = (Gy Di A 
hi + hę i > h; + hi i+2 I zl +2) 
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0<i<n—2, albo w postaci macierzowej 


2 w Ci vı 
us 2 ws c U2 
ug 2 wg cz V3 
g =| ” |, (88) 
Un—2 2 Wn-2 Ch—2 Un—2 
Un1 2 Gl Un—1 
gdzie 
e hi—ı ” hi 
D A EAA AERA 
hi1 + hi hę + hi 


U; = PASZ Tiri). 


Ostatecznie, aby znaleźć współczynniki a;, bi, Ci, d; należy najpierw roz- 
wiązać układ równań (8.8), a potem zastosować wzory (8.7), (8.2), (8.6) 
i (8.5). 

Zauważmy, że macierz układu (8.8) jest trójdiagonalna i ma domi- 
nującą przekątną. Układ można więc rozwiązać kosztem proporcjonal- 
nym do wymiaru n używając algorytmu przeganiania z U. 3.6. Koszt 
znalezienia wszystkich współczynników kubicznej funkcji sklejanej in- 
terpolującej f jest więc też proporcjonalny do n. 


Na końcu oszacujemy jeszcze błąd interpolacji naturalnymi kubicz- 
nymi funkcjami sklejanymi na przedziale |a, b|. Będziemy zakładać, że 
f jest dwa razy różniczkowalna w sposób ciągły. 


Twierdzenie 8.3 Jeśli f € F;y(|a,b]) to 


IF — sę|lcqac) < 5M max (z; = i)”. 


W szczególności, dla podziału równomiernego z; = a + (b — aji/n, 0 < 
i <n, mamy 

b— a2 

IF — ssllCqab) < 5M ( ) . 


n 
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Dowód Wykorzystamy obliczoną wcześniej postać interpolującej funk- 
cji sklejanej sy. Dla z € [x;, 1,41] mamy 


fbl = 7 (2) 
hi 


Me fæ + ( -hlin t20) e-a 


MORETE 

Z rozwinięcia f w szereg Taylora w punkcie z; dostajemy f(x) = f(x;)+ 
Pe — z.) GO S a oraz I aa e eE E 
h.f"(6:)/2. Stąd 


f(z) — s;(z) = f(x) — w(x 


) 
( 


1 1 
- Ue- ap- (LE - (e, 2e) kle) 
—3c;(x zt: xi)? = G1 7 ĉi (x a. xi)’, 
hi 
oraz 
f(x) — sę(z)| < (M+2lej,,| + 6ll)h 
= (M+8 max ([cl)h;. 
1<i<n-1 
Niech teraz max; <j<n_1 |cj| = ||. Z postaci układu (8.8) otrzymujemy 
l] = 2lcj| — |cs|(us + ws) < |uscy q + 20, + WsCs+1| 
UA 
1 
= eseala | E 


a stąd i z (8.9) 
|F(z) — ss(£)| < 5Mh;, 


co kończy dowód. 


Uwagi i uzupełnienia 


U. 8.1 Zamiast terminu funkcje sklejane używa się też często terminów 
splajny (ang. spline-sklejać), albo funkcje gięte. 
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U. 8.2 W Rozdziale 7 pokazaliśmy, że aproksymacja kawałkami wielomia- 
nowa stopnia r jest optymalna (co do rzędu zbieżności) w klasie Fx, (|a,b]), 
wśród wszystkich aproksymacji korzystających jedynie z informacji o warto- 
ściach funkcji w danej liczbie n punktów. Okazuje się, że podobne optymalne 
własności wykazują funkcje sklejane, ale w innych klasach funkcji. 

Niech 


wy(a,b) = { f € W” (a,b): | (FO(a)) dr <M). 


Ustalmy węzły a = 19 < +: < £n = b. Dla f € Wyy(a,b), niech sy będzie 
naturalną funkcją sklejaną interpolującą f w xj, 0 < j < n, a az dowolną 
inną aproksymacją korzystającą jedynie z informacji o wartościach f w tych 
węzłach , tzn. 


az = o(f(zo),..., f(2n)). 


Załóżmy, że błąd aproksymacji mierzymy nie w normie Czebyszewa, ale w 
normie średniokwadratowej, zdefiniowanej jako 


b 
Igllcsch = f Gik 


Wtedy 
sup |f- sflex < sup IIF- asllealab) 
FEW (a,b) fEWṣ; (a,b) 
Aproksymacja naturalnymi funkcjami sklejanymi jest więc optymalna w kla- 
sie Wi; (a,b). 

Można również pokazać, że interpolacja SĘ naturalnymi funkcjami skle- 
janymi na węzłach równoodległych 1; = a + (b — a)j/n, 0 < j < n, jest 
optymalna co do rzędu w klasie Wy,(a,b), wśród wszystkich aproksymacji 
korzystających jedynie z informacji o wartościach funkcji w n+ 1 dowolnych 
punktach, oraz 

* XU eT 
EC IF — sż||ca(ab) x NO”. 
U. 8.3 Tak jak wielomiany, naturalne funkcje sklejane interpolujące dane 
funkcje można reprezentować przez ich współczynniki w różnych bazach. 
Do tego celu można na przykład użyć bazy kanonicznej K,,0<j<n, 
zdefiniowanej równościami 


O iźj, 
kla) = | a 
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przy której s(x) = X, f(c;)K;(x). Baza kanoniczna jest jednak niewy- 
godna w użyciu, bo funkcje K; w ogólności nie zerują się na żadnym pod- 
przedziale, a tym samym manipulowanie nimi jest utrudnione. 

Częściej używa się bazy B-sklejanej Bj, 0 < j < n. W przypadku funkcji 
kubicznych, r = 2, jest ona zdefiniowana przez następujące warunki: 


B;(x;) — 1, dla 0 < J < n, 
B;(x) = 0, dla x < £zj-2,j > 2, oraz dla z > £j42,j < n -— 2. 


Dla Bo i Bı dodatkowo żądamy, aby 
By (zo) = 0 = By(zo), Bı(zxo) = 0, 
a dla By_1 i Bn podobnie 
B! _ (£n) = 0 = B! (z£n), Bn- (tn-1) =0; 


Wtedy B; nie zeruje się tylko na przedziale (1;_2, xj+2), Wyznaczenie współ- 
czynników rozwinięcia w bazie 4B;);_9 funkcji sklejanej interpolującej f wy- 
maga rozwiązania układu liniowego z macierzą trójdiagonalną (B;(z;))7;—0> 
a więc koszt obliczenia tych współczynników jest proporcjonalny do n. 


U. 8.4 Oprócz naturalnych funkcji sklejanych często rozpatruje się też okre- 
sowe funkcje sklejane. Są to funkcje 5 : R > R spełniające warunki (i), (ii) 
z Rozdziału 8.1, oraz warunek: 


(iii)? 5Ó) jest dla 0 < i < r — 1 funkcją okresową o okresie (b — a), tzn. 
56(x) = 5® (x + (b — a)), Yz. 


Klasę okresowych funkcji sklejanych rzędu r oznaczymy przez Ś,. Funkcje 
te mają podobne własności jak naturalne funkcje sklejane. Dokładniej, niech 


W'(a,b) = ff E W'(a,b) : fÓ(a)=f(b), O<i<r-1), (8.9) 
tzn. W” (a,b) jest klasą funkcji z W” (a, b), które można przedłużyć do funkcji, 
krórych wszystkie pochodne do rzędu r — 1 włącznie są (b — a)-okresowe na 


R. Wtedy dla dowolnej funkcji f € W” (a,b) zerującej się w węzłach zj, oraz 
dla dowolnej $ € 5, mamy 


M * Oaz) dz = 0. (8.10) 
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Jest to odpowiednik Lematu 8.1 w przypadku okresowym. W szczególno- 
ści wynika z niego jednoznaczność rozwiązania zadania interpolacyjnego dla 
okresowych funkcji f (tzn. takich, że f(a) = f(b)), jak również odpowiednia 
własność minimalizacyjna okresowych funkcji sklejanych. Dokładniej, jeśli 
f € W" (a,b) oraz Sp E Š, interpoluje f w węzłach Tj, 0<J <n, to 


T (FO) de > a (EP (a)) da. 


U. 8.5 Klasyczne zadanie aproksymacyjne w przestrzeniach funkcji defi- 
niuje się w następujący sposób. 

Niech F będzie pewną przestrzenią liniową funkcji f : [a,b] > R, w 
której określona została norma | : ||. Niech V, C F będzie podprzestrzenią 
w F wymiaru n. Dla danej f € F, należy znaleźć funkcję vę € F taką, że 


If- vgl = min I£ — vll- 


Okazuje się, że tak postawione zadanie ma rozwiązanie v;, choć nie zawsze 
jest ono wyznaczone jednoznacznie, zob. Ów. 8.6. 


Jako przykład, rozpatrzmy F = W” (a,b). Utożsamiając funkcje fi, f2 € 
W”'(a,b) takie, że f4(z) — falx) € II,_;, zdefiniujemy w W”(a,b) normę 


I/I = T [| OOP a 


Dla ustalonych węzłów a = 19 < +- < £n = b, niech 


Vn+1 = Sr 


będzie podprzestrzenią w W” (a,b) naturalnych funkcji sklejanych rzędu r 
opartych węzłach zj, 0 < j < n. Oczywiście dims, = n + 1, co wynika z 
jednoznaczności rozwiązania w 6, zadania interpolacji. Okazuje się, że wtedy 
optymalną dla f e W'(a,b) jest naturalna funkcja sklejana sy interpolująca 
f w węzłach zj, tzn. 


If- szl = min If -— sll. 


Rzeczywiście, ponieważ norma w przestrzeni W” (a,b) generowana jest 
przez iloczyn skalarny 


b 
TESE f (a) (E) dz, 
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jest to przestrzeń unitarna. Znane twierdzenie mówi, że w przestrzeni uni- 
tarnej najbliższą danej f funkcją w dowolnej domkniętej podprzestrzeni V 
jest rzut prostopadły f na V, albo równoważnie, taka funkcja vę € Vn+1, że 
iloczyn skalarny 

d=upiU = 0, Ww eV. 


W naszym przypadku, ostatnia równość jest równoważna 
b 
l (f - vs) (x)(x) dx = 0, Vs ES, 
a 


To zaś jest na mocy Lematu 8.1 prawdą gdy v; interpoluje f w punktach 
Tj, czyli vf = Sy. 

Dodajmy jeszcze, że nie zawsze interpolacja daje najlepszą aproksymację 
w sensie klasycznym, zob. Cw. 8.7. 
Cwiczenia 
Ćw. 8.1 Dla z€ R, niech 


z z>0, 
z+ = max(0, z) = i | W) 


Pokazać, że każdą funkcję sklejaną s rzędu r można przedstawić w postaci 
2r—1 
s(x) = w(x) + X a;((a — ts) e 


gdzie a; są pewnymi współczynnikami rzeczywistymi, a w € Llo,_1. Ponadto, 
jeśli s jest naturalna, to w € II,_1, a współczynniki a; spełniają równania 


n 
X ajzi = 0, 0<4<r-1. 
j=0 


Ćw. 8.2 Niech h > 0 i c e R. Wyznaczyć współczynniki kubicznej funkcji 
sklejanej s opartej na pięciu węzłach —2h, —h, 0, h, 2h i spełniającej dodat- 
kowo następujące warunki interpolacyjne: 


s(0) = c, s™®(+2h) = 0,  k=0,1,2, 


102 ROZDZIAŁ 8. INTERPOLACJA FUNKCJAMI SKLEJANYMI 


Ćw. 8.3 Wykorzystując Ćw. 8.2 wykazać, że dla dowolnego układu węzłów 
aŚaq <:''xa, Śb mamy 


Myb—a)2 
n ; 1 A= 1 < >> Z ź 
maxq ||f'||cqax) : f € Fu(la,b]), f(x) 0,1 <i<n) > w ga) 


Stąd i z U. 7.6 wywnioskować, że dla każdej aproksymacji a; korzystającej 
jedynie z wartości f w n punktach, 


FREE Z 
max —a — : 

jeFy(abp RRE AARET 

Ćw. 8.4 Pokazać, że funkcje Bj, 0 < j < n, zdefiniowane w U. 8.3 istnieją 
i tworzą bazę w przestrzeni S,. 


Ćw. 8.5 Pokazać jednoznaczność rozwiązania zadania interpolacyjnego w 
przypadku okresowych funkcji sklejanych (zob. U. 8.4), oraz własność mini- 
malizacyjną (8.9). 

Wskazówka. Zastosować technikę dowodową podobną do tej z przypadku 
naturalnych funkcji sklejanych. 


Ćw. 8.6 Pokazać, że ogólnie postawione zadanie aproksymacyjne z U. 8.5 
ma rozwiązanie. 

Wskazówka. Zauważyć, że jeśli v istnieje to ||v;|| < 2||f||, oraz wykorzystać 
ciągłość odwzorowania v > ||f —v||, v € V, oraz fakt, że dimV < oo. 


Ćw. 8.7 Niech F = C([a,b]) z normą Czebyszewa, albo F = C>(a,b) z 
normą średniokwadratową. Niech V,+1 = In. Uzasadnić na przykładach, że 
w ogólności nie istnieją węzły xj, 0 < j < n, o następującej własności: dla 
każdej funkcji f € F, wielomian wę € II, interpolujący f w węzłach z; jest 
elementem optymalnym w sensie aproksymacji klasycznej, zdefiniowanej w 
U. 8.5. 


CE 


Rozdział 9 


Całkowanie numeryczne 


Zajmiemy się teraz zadaniem całkowania numerycznego. Polega ono na 
obliczeniu (a raczej przybliżeniu) całki oznaczonej 


gdzie —0o < a < b < +œ, a f należy do pewnej klasy F funkcji 
rzeczywistych określonych i całkowalnych w sensie Riemanna na całym 
przedziale [a,b]. 

Będziemy zakładać, że mamy możliwość obliczania wartości funkcji 
f, a w niektórych przypadkach również jej pochodnych, o ile istnieją. 
Dokładna całka S(f) będzie więc w ogólności przybliżana wartością 
A(f), która zależy tylko od wartości f i ew. jej pochodnych w skończo- 
nej liczbie punktów. 


9.1 Co to są kwadratury? 


Kwadraturami nazywamy funkcjonały liniowe Q : F — R postaci 


n 


QC) = X afla), 


i=0 


albo ogólniej 
k Ni —1 


AS S a a; (9.1) 


i=0 j=0 
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gdzie x, są punktami z [a,b], a a, (albo a,,) są pewnymi współczyn- 
nikami rzeczywistymi. Zauważmy, że obliczenia kwadratur są dopusz- 
czalne w naszym modelu obliczeniowym, mogą więc służyć jako sposób 
przybliżania całki. 

Jeden z możliwych sposobów konstrukcji kwadratur jest następu- 
jący. Najpierw wybieramy węzły x; (pojedyncze lub wielokrotne), bu- 
dujemy wielomian interpolacyjny odpowiadający tym węzłom, a na- 
stępnie całkujemy go. Ponieważ postać wielomianu interpolacyjnego 
zależy tylko od danej informacji o f, otrzymana w ten sposób wartość 
też będzie zależeć tylko od tej informacji, a w konsekwencji funkcjonał 
wynikowy będzie postaci (9.1). Są to tzw. kwadratury interpolacyjne. 


Definicja 9.1 Kwadraturę Q! opartą na węzłach o łącznej krotności 
n + 1 nazywamy interpolacyjną, jeśli 


QG) = f wod, 


gdzie w; jest wielomianem interpolacyjnym funkcji f stopnia co najwy- 
żej n, opartym na tych węzłach. 


Współczynniki kwadratur interpolacyjnych można łatwo wyliczyć. 
Rozpatrzmy dla uproszczenia przypadek, gdy węzły są jednokrotne. 
Zapisując wielomian interpolacyjny w postaci jego rozwinięcia w bazie 
kanonicznej Lagrange'a l, (zob. (6.3)), otrzymujemy 


= e aas > Fa) [ ula) da. 


a stąd i z postaci li, 


z = (© — zo): (£ — Ti—1)(£ — zj) (£ — Tn) A 
zB £o) o. zj Ti—1)(qi E Z) EA (zi ej Tn) | 
0<i<n. 
Podamy teraz kilka przykładów. 
Kwadratura prostokątów jest oparta na jednym węźle zo = (a + b)/2, 


QU) = b= a)f( 


+b 
— ). 
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Kwadratura trapezów jest oparta na jednokrotnych węzłach zo = a, 
zı = b i jest równa polu odpowiedniego trapezu, 


AP = TA) = = (Fla) + F6)). 


Kwadratura parabol (Simpsona) jest oparta na jednokrotnych węzłach 
£o = a, 1, = b, £2 = (a +b)/2, i jest równa polu pod parabolą interpo- 
lującą f w tych węzłach, 


b—a 


(ro +aj (=) +70). 


Zauważmy, że kwadratury trapezów i parabol są oparte na węzłach 
jednokrotnych i równoodległych, przy czym zo = ai £n = b. Ogólnie, 
kwadratury interpolacyjne oparte na węzłach równoodległych z, = a + 
(b= aji/n,0O<i<n, nazywamy kwadraturami Newtona-Cotesa. 


Qz(f) = P(F) = 


9.2 Błąd kwadratur interpolacyjnych 


Zajmiemy się teraz błędem kwadratur interpolacyjnych. Przypomnijmy, 
że Fyx,(|a,b|) oznacza klasę funkcji (n + 1) razy różniczkowalnych w 
sposób ciągły i takich, że |f"*V(x)| < M, Va. 


Twierdzenie 9.1 Niech Q? będzie kwadraturą interpolacyjną opartą 
na (jednokrotnych lub wielokrotnych) węzłach x, 0 < i < n. Jeśli f € 
Fy(la,b]) to 


zai g)>. (9.2) 


W klasie F% (Ja, b|), maksymalny błąd kwadratury Q! wynosi 
M b 
I = 
p ISA- = mog J, (E706-0)-—(6-2,)|dz. 


FEF, (lab) (n 


Dowód Korzystając ze znanego nam już wzoru na błąd interpolacji 
wielomianowej z Lematu 7.1, mamy 


SP) - QH) = f(E- a)(e- 21) (8 —2n)f(@o 81+- Em 2) de. 
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Stąd, jeśli f € Fxy,((a,b|) to 


M 
(n+ 1)! 


SF) - Q (Pl < [e-a dz = (b — a)" *? 


M 
a (n+ 1)! 


Ograniczenie górne w dokładnej formule na błąd w klasie Fx,([a, b]) 
wynika bezpośrednio z oszacowania (9.2). Aby pokazać ograniczenie 


dolne zauważmy, że dla funkcji g takiej, że g(”*") przyjmuje na prze- 
działach (a, £o), (£o, £1), ..., (£n, b0) naprzemiennie wartości M i -M 
mamy 


M _p 
5) - 0] = zyj, |E-206-0)--(G6-z,|dz. 


Co prawda, g nie jest w Fy|a,b|, ale może być dla dowolnego e > 0 
przybliżana funkcjami f; € F% ([a,b]) w ten sposób, że całka 


F l(e-m)--c-2(F- 0) a) de < e(n+ 1). 


Zapisując fe = g + (fe — g) mamy 
IS(f.) — Q4| < IS(9) — Q'(9)| a — 9) — Q'(F:— 9)| 
< EA CZE )---(1 — 2n)|dr + e 
T (n+1]!Ja 0 É í 


co wobec dowolności € daje dowód twierdzenia. 


W szczególnych przypadkach kwadratur trapezów T i parabol P 
możemy otrzymać innego rodzaju formuły na błąd. 


Twierdzenie 9.2 (i) Jeśli f € C”([a,b]) to dla kwadratury trapezów 
mamy 


— 3 
sA) - Td) E 
(ii) Jeśli f e C(f[a,b]) to dla kwadratury parabol mamy 
S 0) pa 
SA) - PI) = — zy "(63 


(M 52 < la, b)). 
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Dowód (i) Ze wzoru (9.2) 


s) -TO= [ 6-d6G-dilabr)dz. 


Ponieważ funkcja z > f(a,b, x) jest ciągła, a wielomian (x — a)(x — b) 
przyjmuje jedynie wartości nieujemne, można zastosować twierdzenie o 
wartości średniej dla całki, aby otrzymać 

b 


SM - TA) = Fabo) | (c= a(e= b) de 


a 


136) 6— a)? 
2! 6 


dla pewnych c, 4 € [a,b]. 

(ii) Niech wf € Ib i wz E€ Ilg będą wielomianami interpolacyjnymi 
funkcji f odpowiednio dla węzłów a, b, (a + b)/2 oraz a, b, (a +b)/2, (a + 
b)/2. Wtedy 


wę3(1) = węo(1) + i (ab, E, e- ay(a- Ehei), 
Wobec ; leż 

| E-a(= ) (a — b) dz =0 
mamy 


P() = f wjafa)dr = [ wata) de. 


Stąd i ze wzoru na błąd interpolacji Hermite'a otrzymujemy 
b 


SD - PO) = | (-wa)la)da 
b a + by2 a+b a+b 
= i (z -a)(x — ) (z — b)f (a,b, yoa o dz. 
Ponieważ wielomian (x — a)(x — (a + b)/2)°(x — b) jest niedodatni na 
la, b], możemy znów zastosować twierdzenie o wartości średniej. Mamy 


b b 
S(f) — P(F) = (ab = — „e) 


[c-a(- O) (a — 1) dx 
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co kończy dowód. 


9.3 Kwadratury złożone 


Podobnie jak w przypadku zadania interpolacji chcielibyśmy, aby błąd 
kwadratur malał do zera, gdy liczba węzłów rośnie do nieskończoności. 
Można to osiągnąć stosując np. kwadratury złożone. Są to kwadratury, 
które powstają przez scałkowanie funkcji kawałkami wielomianowej in- 
terpolującej f. 

Prostym przykładem kwadratury złożonej jest suma Riemanna, 


n 


PEOS T t)f(2:), 


i=0 


gdzie a = tę < ti <- < ty = b oraz x; € [ti ti+1]- Jeśli średnica 
podziału, maxo<i<n(ti — ti-1), maleje do zera to limm_o Q(f) = S(f). 
Będziemy rozpatrywać kwadratury złożone postaci 


GH) = [| wy(a)dz, 
gdzie wy jest kawałkami wielomianem z Rozdziału 7.3. To znaczy, dla 
danego n kładziemy t; = a + (b — aji/k, 0 < i < k, a następnie dla 
każdego i wybieramy dowolne węzły £i; € [ti-1, ti], 0 < j < r. Wtedy 
w; jest na każdym przedziale wielomianem interpolacyjnym funkcji f 
stopnia co najwyżej r opartym na węzłach zx; ;. Kwadratura Q korzysta 
z węzłów o łącznej krotności n < k(r + 1). 


Twierdzenie 9.3 Błąd kwadratury złożonej Q(f) w klasie F;y([a,b]) 
jest ograniczony przez 


E (b — a)"*? M ląr+1 
ib POOS A= NSG2 


feFy,([a,b]) 
gdzie 
M(r dE me = ajr 


AE (r + 1)! 
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Dowód Twierdzenie to jest bezpośrednim wnioskiem z Twierdzenia 
9.1. Mamy bowiem 


s0-QAD < Ef 0- wydldz 
k ,b=ar+2 M (b-a? M 
S >| ( a , 


r+1) kt (r+1)" 


co kończy dowód. 


W klasie F;,([a,b]), błąd kwadratur złożonych jest rzędu n""*V, 


czyli jest taki sam jak błąd interpolacji kawałkami wielomianowej. I 
tak jak przy interpolacji można pokazać, że błąd każdej innej metody 
całkowania korzystającej jedynie z wartości funkcji w n punktach nie 
może w klasie F;,([a,b]) maleć szybciej niż n""*, zob. U. 9.1. Podane 
kwadratury złożone mają więc optymalny rząd zbieżności. 

Zajmiemy się teraz błędem szczególnych kwadratur złożonych, mia- 
nowicie złożonych kwadratur trapezów Ty i parabol P,. Powstają one 
przez zastosowanie na każdym przedziale |£,_4, t;| odpowiednio kwadra- 
tur trapezów T i parabol P. Jak łatwo się przekonać, 


4 = 8 (OLE | LD). 


Pf) = yi | EÈ) + ŻY). 
Twierdzenie 9.4 (i) Jeśli f € C”([a,b)) to 
SA) - TA) = ze Og), 
(ii) Jeśli f € C®([a,b]) to 
SA) - BG) = E 
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Dowód Dla kwadratury trapezów mamy 


SU) - RH) = -D E (0) 
„be iti oya Oai 
= -CS Oo) = -E Og), 


i=l 


a dla kwadratury parabol podobnie 


SH) - BA) = -X zza (A) 
_(b= a)? k e. (b — a)” „qą 


2280k! >. k E= soa (6). U 


Kwadratura parabol ma więc optymalny rząd zbieżności nie tylko 
w klasie F? ([a,b]), ale też w Fẹ (la, b]). 


9.4 Przyspieszanie zbieżności kwadratur 


W praktyce często stosuje się obliczanie kwadratur poprzez zagęszcza- 
nie podziału przedziału [a,b]. Na przykład, dla złożonej kwadratury 
trapezów zachodzi następujący wygodny wzór rekurencyjny: 


Ty = ; (0 + a VË) (9.3) 


Pozwala on obliczyć Tx,(f) na podstawie T,(f) poprzez “doliczenie” 
wartości funkcji w punktach *gęstszej” siatki. W ten sposób możemy 
obserwować zachowanie się kolejnych przybliżeń 1>.(f) (s > 0) całki 
S(f). Jest to szczególnie istotne wtedy, gdy nie mamy żadnej informa- 
cji a priori o ||f"||c(favj); a przez to nie potrafimy oszacować liczby n 
węzłów, dla której osiągniemy pożądaną dokładność, zob. U. 9.2. 

Jeśli funkcja jest więcej niż dwa razy różniczkowalna to użycie złożo- 


nych kwadratur trapezów zdaje się tracić sens. Wtedy istnieją przecież 
kwadratury, których błąd maleje do zera szybciej niż n**. Okazuje się 
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jednak, że kwadratury Tą mogą być podstawą dla prostej rekurencyjnej 
konstrukcji innych kwadratur posiadających już optymalną zbieżność. 
Konstrukcja ta bazuje na następującym ważnym lemacie. 


Lemat 9.1 (Formuła Fulera-Maclaurina) E 
Dla funkcji f e CC™+2([a,b]), błąd złożonej kwadratury trapezów Th, 
wyraża się wzorem 


SA) - DAF) = E ah” (fe) — fQ-V(a)) 


i=1 
+ Cmpa h? t? (b — a) FOTO Enk), 
gdzie h = (b — a)/k, Em € [a,b], a ci są pewnymi stałymi liczbowymi. 
Mamy cı = —1/12, co = —1/720 i, ogólnie, c; = B;/(2i)!, gdzie B; są 
tzw. liczbami Bernoulliego. O 


Dowód tego lematu pominiemy. 


Formułę Eulera-Maclaurina można przepisać w postaci 
sA - D) = LAK + aal) eTA, 
i=1 


gdzie c(?(f) = a(b — až (fD) — fE-D(a)), 1 < i < m, oraz 
chik(f) = Cmq1(b — a) fOD Enyi). Zauważmy przy tym, że 
jeśli f e Fy*'([a,b]) to współczynniki e ddd0) są wspólnie ograni- 
czone przez Cm_1(b — a) ™ t? M. 

Definiując teraz kwadraturę 


dla f e C® (fa, b]) mamy 


DER E ASCE) — Tal H — (S(F) TC) 


3 
aa Aa E a 
~ 3 \ 4k? ` 42h 3\ kK ` kê 
20) 


po 


112 ROZDZIAŁ 9. CAŁKOWANIE NUMERYCZNE 


gdzie E) = (ide) — (1/3)c$) (7) i jest wspólnie ograniczone 
dla f € F} (fa, b]). Kwadratura T} ma więc optymalny w F} (fa, b]) rząd 
zbieżności k74. Proces ten można kontynuować dalej tworząc kolejne 
kwadratury o coraz to wyższym rzędzie zbieżności. Dokładniej, połóżmy 
TY(f) = T,( f) oraz, dla s > 1, 


srns—l __ rns—1 


(9.4) 


Wtedy, dla f € Fy"*'([a,b]), rząd zbieżności kwadratury T} wynosi 
k-©m+2), Rzeczywiście, sprawdziliśmy, że jest to prawdą dla m = 0, 1. 
Niech m > 2. Postępując indukcyjnie ze względu na s = 1,2,...,m 
mamy 


Sj) — ej) = FEG) — Te M (SA Sn UD 


J (e (> a O(P + daue ea) 


a s— —24 s— —(2m 1 
z p PAPET + czę e(P)ET* o RE 
=.) 02006 T ŚR, 
1=s+1 


ponieważ współczynniki przy k7? redukują się. c(*(f) są tutaj pew- 


nymi nowymi stałymi, a A „(f) może być w klasie F?™+! (Ja, b]) ogra- 
niczona przez stałą niezależną od f. Ostatecznie, dla s = m mamy więc 


SDE =d „(s 
i w klasie Fyy”*'([a,b]) 
BORMAAIESZT ANA 


dla pewnej stałej c, niezależnej od f. 

Zauważmy jeszcze, że I)? wykorzystuje n = k2™ + 1 wartości f 
w punktach równoodległych na [a,b] co oznacza, że w terminach n 
rząd zbieżności wynosi też n"©”"*), a więc jest optymalny w klasie 
Fir (la, b). 
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Kwadratury Tf nazywane są kwadraturami Romberga. Dla danej 
funkcji f, można je łatwo konstruować budując następującą tablicę trój- 
kątną: 


rf) 
nA) TO | 

TD DO RA m 
A TO RO BO) (2) 


RO TL Ta) Bo —: FAD, 


w której pierwsza kolumna jest tworzona indukcyjnie zgodnie ze wzorem 
(9.3), a kolejne zgodnie z (9.4). 


Uwagi i uzupełnienia 


U.9.1 Pokażemy, że jeśli dowolny algorytm A daje przybliżoną wartość 
całki S(f) = f A f(x)dx na podstawie jedynie wartości lub pochodnych funk- 
cji f w n różnych punktach, tzn. 


Aea "laka sa (6) 
dla pewnych z; € [a,b], to błąd takiego przybliżenia w klasie Fyy([a,b]), 


sup  |S(f) — A(f)| GM n n, (9.6) 
jEFy, (lad) 


gdzie C4 jest pewną stałą dodatnią niezależną od A i n. 
W tym celu, podobnie jak w U. 7.2, wybierzmy dowolną nieujemną funk- 
cję Y : R — R spełniającą warunki: 


1 y € FI (R), 
2. p(x) = 0 dla x æ [0,1]. 


Oznaczmy zo = 0, £n+1 = b, hi = Gy — zi, 0 < i < n, oraz zdefiniujmy 


funkcje 


dia) = MK A) 
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Wtedy dla g,(x1) = —go(x) = > 5-9 %;(1) mamy gı E€ Fy(la,b|) oraz infor- 
macja o gı jest zerowa, gP (z5) = 0 da0 <p<r+1,1<j< n. Stąd 
o = $(0,...,0) jest aproksymacją dla całek obu funkcji w i 

—— 


n(r+2) 
maxt|S(a) — dol, |S(92) — gol} 2 |S(m) — S(92)1/2 = Slon) 
z | > s" Roa | (20 dz 


i=0)* Ti-1 = i—1 


1 n 
u f y(a) dr + Y KT? 
0 i=0 


Zauważmy jeszcze, że suma $>% o hy” jest minimalna dla h; = (b—a)/(n+1), 
Vi, a stąd 


— a)" +2 
sim) > MEIG f vta)dz. 


Nierówność (9.6) zachodzi więc z C1 = (2(b — a))"F! Ją y(x)dz. 

U. 9.2 Jeśli f e CÓ([a,b)) i f'(a) £ f'(b) to dla dużych k mamy 
S(f) — T(f)  clf)k*, 

gdzie c(f) = —(b — a)? (f'(b) — f'(a))/12. Stąd 

TŁ(F) — Thal) = (SCF) — Tkja(f)) — (SUP) — Thl) 


) _ 3c0) B- 
k2 (k/2)2 k 


a więc 
SW) - DU) = z Gal) — Gl1)). 


gdy k — oo. Ostatnia przybliżona równość może sugerować kryterium koń- 
czenia obliczeń, gdy chcemy otrzymać przybliżenie całki S(f) z zadaną do- 
kładnością e > 0. Mianowicie, obliczając kolejne kwadratury T>:(f) dla 
s = 0,1,..., sprawdzamy jednocześnie, czy 


[5—(f) — T>s(f)| < 3e. 


Jeśli ten warunek jest spełniony dla kilku kolejnych s to z dużym prawdo- 
podobieństwem możemy stwierdzić, że |S(f) — T>s(f)| < e. 

Oczywiście, podobne kryteria kończenia obliczeń mogą być konstruowane 
dla innych ciągów kwadratur, jeśli tylko znamy teoretyczne zachowanie się 
błędu. 
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Cwiczenia 


Ćw. 9.1 Pokazać, że w kwadraturze Newtona-Cotesa opartej na n + 1 rów- 
noodległych węzłach z; = a + (b — aji/n, współczynniki a; dane są wzorami 


da; = 


GE | 28-0--6-6-1e- (+0) (en) da: 


ni!(n — i)! 


Ćw. 9.2 Pokazać, że jeśli f e C)([a,b)) to dla kwadratury prostokątów 
Qo(F) = F((a + b)/2)(b — a)/2 mamy 
(b— a)? 


S) — Golf) = —ą tP 6o), 


(ĉo € [a,b]), a w konsekwencji 


AŻ — Qo(f)| = =. 


Błąd w klasie F} ([a,b]) jest więc dla kwadratury prostokątów dwa razy 
mniejszy niż dla kwadratury trapezów. Wywnioskować stąd, że złożona kwa- 
dratura prostokątów, 


ma optymalny rząd zbieżności nie tylko w klasie F$y([a, b]), ale też w klasie 
Fy(la,b]). 


Ćw. 9.3 Rozpatrzmy kwadratury interpolacyjne oparte na dwóch węzłach 
zo, £1 E [a,b]. Pokazać, że wśród tych kwadratur najmniejszy błąd w klasie 
F+([a,b]) jest osiągany przez kwadraturę 


Q!) = = ËH | (7). 


a jej błąd 


b— 3 
sup [SP - QPI = Z. 
fEFY ([a,b]) 
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Ćw. 9.4 Uzasadnić, że jeśli kwadratura interpolacyjna Q! oparta jest na 
węzłach Czebyszewa 


a+b b-a 2i +1 ; 
e <i< 
Ti 9 F 9 cos (737): aar, 
to dla f € Fy([a,b]) błąd 
M(b — a)"t? 


SA - GD] < T 


a dla odpowiadającej węzłom Czebyszewa kwadratury złożonej 


IS(A) — BAI < 


M (b — a)"** ok 
22r+1(r + 1)! Nk 


Ćw. 9.5 Rozpatrzmy kwadraturę 


= (b — a)? 


GA) = KA) -S O) - FA) 


Pokazać, nie korzytając z formuły Eulera-Maclaurina, że jeśli f € F;,([a,b]) 
to błąd |S(f) — Qz(f)| zbiega do zera co najmniej tak szybko jak k”*, gdy 
k — oo. W szczególności, dla funkcji spełniających f'(a) = f'(b) rząd zbież- 
ności złożonej kwadratury trapezów wynosi k”*. 


Ćw. 9.6 Opracować ekonomiczny algorytm obliczania ciągu kolejnych zło- 
żonych kwadratur parabol P>s(f) dla s=0,1,2,.... 


Ćw. 9.7 Opracować kryterium kończenia obliczeń analogiczne do tego z U. 
9.2, dla ciągu złożonych kwadratur parabol Pəs(f), s=0,1,.... 


Ćw. 9.8 Pokazać, że drugą kolumnę tabeli (9.5) kwadratur Romberga two- 
rzą złożone kwadratury parabol, tzn. 


Alf) = PAT 


Ćw. 9.9 Opracować program obliczający wartość T;(f) kwadratury Rom- 
berga, korzystając ze wzorów (9.3) i budując tabelę (9.4). 


CE 


Rozdział 10 


Całkowanie a aproksymacja 


W poprzednim rozdziale badaliśmy kwadratury ze względu na kryte- 
rium błędu w pewnych klasach funkcji. W tym rozdziale zajmiemy się 
badaniem kwadratur ze względu na inne, bardziej klasyczne kryterium, 
którym jest rząd kwadratury. Będziemy przy tym rozpatrywać nieco 
ogólniejsze zadanie całkowania, tzw. całkowanie z wagą, 


gdzie waga p jest prawie wszędzie dodatnia i jest całkowalna, 


b 
| |plr)|dr < +. 


Dopuścimy również możliwość nieskończonych przedziałów całkowania, 
a więc oo <a < b < +00. 


10.1 Rząd kwadratury 


Zaczniemy od definicji. 


Definicja 10.1 Rzędem kwadratury Q nazywamy taką liczbę r, że 


(i) kwadratura Q daje dokładną wartość całki dla wszystkich wielomia- 
nów stopnia mniejszego niż r, 


Q(w) = S,(w), Yw € Il, 


117 
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(i) istnieje wielomian stopnia r dla którego kwadratura nie jest równa 
całce, 


+w* € Il, Q(w*) £ S(w*). 


Rząd kwadratury Q oznaczymy przez rz((0). Dla przykładu, za- 
łóżmy, że p = 1 i przedział [a,b] jest skończony. Wtedy rząd kwa- 
dratury prostokątów wynosi 2, tak jak rząd kwadratury trapezów T, 
natomiast rząd kwadratury parabol P wynosi 4. Rzeczywiście, to wy- 
nika bezpośrednio ze wzorów na błędy tych kwadratur danych w Ćw. 
9.2 i Twierdzeniu 9.2. Ogólniej, mamy następujące twierdzenie, które 
przytaczamy bez dowodu. 


Twierdzenie 10.1 Niech waga p = 1 i niech przedział [a,b] będzie 
skończony. Rząd kwadratury Newtona-Cotesa Q)*, tzn. opartej na jed- 
nokrotnych węzłach równoodległych x; = a + (b — aji/n,0<i<n, 
wynosi 


SER]. BI dla n nieparzystych, 
rz(Q ) = n+2 dla n parzystych. 


Jasne jest, że interesują nas kwadratury o jak najwyższym rzędzie. 
Chcielibyśmy wiedzieć, jaki jest maksymalny rząd kwadratury korzy- 
stającej z ustalonej liczby węzłów i jak konstruować kwadratury o mak- 
symalnym rzędzie. Jak przekonamy się później, nie jest to tylko aka- 
demickie pytanie, bowiem kwadratury o maksymalnym rzędzie mają 
również dobre własności ze względu na błąd. 


Częściową odpowiedź na pytanie o maksymalny rząd kwadratur 
daje następujący lemat. 


Lemat 10.1 Niech Q będzie kwadraturą opartą na węzłach o łącznej 
krotności n + 1. 

(i) Jeśli rz(Q) > n+1 to Q jest kwadraturą interpolacyjną. 

(ii) Jeśli Q jest kwadraturą interpolacyjną to 


n+1 <rz(Q) <2(n+1). 
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Dowód (i) Dla dowolnej funkcji f, niech wę będzie wielomianem inter- 
polacyjnym dla f stopnia co najwyżej n, opartym na tej samej informa- 
cji o f, z jakiej korzysta kwadratura Q. Wtedy mamy Q(f) = Q(w;), 
a ponieważ rz(ęQ) > n + 1 to również Q(w;) = S,(wyę). Stąd 


Q) = Sp(wy), 


a to oznacza, że Q jest interpolacyjna. 
(ii) Jeśli Q jest interpolacyjna to dla każdego w € II, mamy Q(w) = 
S (w), bo w jest dla siebie wielomianem interpolacyjnym. Stąd rz(Q) > 
n+l. 

Dla dowodu drugiej nierówności załóżmy, że Q korzysta z węzłów 
To, Z1,- --, Zn (być może powtarzających się). Wtedy dla wielomianu 


w,(z) = (z — zo) (z — 24)” *--(z— 2)” 
mamy w; € Ilo, +1) Oraz 5,(w,) > 0, bo wą jest prawie wszędzie do- 
datni. Z drugiej strony, Q(w,) = 0, bo informacja o wy jest zerowa. 
Stąd 

Q(wi) # Sow) 


i w konsekwencji rz(Q) < 2(n +1). © 


Lemat 10.1 mówi nam, że w poszukiwaniu kwadratur o najwyż- 
szym rzędzie możemy ograniczyć się do kwadratur interpolacyjnych. 
Ponadto, maksymalny rząd jest na pewno nie mniejszy niż n + 1 i na 
pewno nie większy niż 2(n + 1). Aby jednak odpowiedzieć na pytanie 
jaki jest rzeczywisty maksymalny rząd, musimy odwołać się do pewnych 
faktów z teorii aproksymacji w przestrzeniach unitarnych, a ściślej, do 
własności ciągów wielomianów ortogonalnych. 


10.2 Ciągi wielomianów ortogonalnych 


Waga p definiuje nam przestrzeń £2 (a,b) funkcji ciaągłych o warto- 
ściach rzeczywistych, określonych na przedziale (a, b) i takich, że całka 
J? £?(«)p(a) dz istnieje i jest skończona. W £2 (a, b) wprowadzimy ilo- 
czyn skalarny 


(6,9) = f FOE) de 
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i generowaną przez ten iloczyn normę 


IAI = VED = y f Pæ) de 


W ten sposób, L2 (a,b) jest przestrzenią unitarną i można mówić o 
ortogonalności (prostopadłości) elementów w tej przestrzeni. 


Definicja 10.2 Ciąg wielomianów pp, k = 0,1,2,..., nazywamy ciq- 
giem wielomianów ortogonalnych w przestrzeni La (a,b) jeśli 


(i) wielomian py jest stopnia dokładnie k, 


deg pe = k,  Vk>0, 


(ii) wielomiany te są wzajemnie ortogonalne w przestrzeni Lə (a,b), 


tzn. 
b 


(pup) = | plop(ap(o)dr=0, ii 


Pokażemy, że ta definicja ma sens. 


Lemat 10.2 Ciąg wielomianów ortogonalnych 4pkjk>o w przestrzeni 
unitarnej Lo pla, b) istnieje. Jeśli dodatkowo założymy, że py(z) = x" + 
..., tzn. współczynnik przy najwyższej potędzie wielomianu pę wynosi 
1, to ciąg taki jest wyznaczony jednoznacznie. 


Dowód Wielomiany py można skonstruować stosując proces ortogo- 
nalizacji Grama-Schmidta, tzn. ortogonalizując dowolny ciąg wielomia- 
nów (w; |x>o spełniający deg wy = k. Rzeczywiście, weźmy np. wy(z) = 
x”. Jeśli położymy po = wy i dalej indukcyjnie 


j=0 (Pr Ba) 


Pj, k > T; 


to, jak łatwo bezpośrednio sprawdzić, (pz }x>0 będzie ciągiem wielomia- 
nów ortogonalnych i p(x) = zë +.... 

Aby pokazać jednoznaczność załóżmy, że istnieje inny ciąg {qk}k>0 
wielomianów ortogonalnych, w którym współczynniki przy najwyższych 
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potęgach wynoszą 1. Ponieważ deg py = k, Vk, wielomiany po, p1, ... „Pk 
tworzą bazę ortogonalną w Il; (ze względu na iloczyn skalarny (-,:)), a 
stąd 


Z drugiej strony, dla 0 < j < k—1, mamy (qr, pj} = 0, bo ortogonalność 

q do wielomianów bazowych qj, 0 < j < k — 1, jest równoważna 

ortogonalności do całej przestrzeni IIķ—1, a w szczególności do p;. Stąd 
(dk, Pr) 


k= Pk, 
(Dk, Pk) 


a ponieważ współczynniki przy x* w wielomianach q i py są takie same, 
to qk = pk. To kończy dowód. 


Pokazaliśmy, że ciąg wielomianów ortogonalnych jest wyznaczony 
jednoznacznie z dokładnością do współczynników przy najwyższych po- 
tęgach z. Podamy teraz kilka przykładów takich ciągów. 

Wielomiany Legendre a Ly. 

Jest to najbardziej naturalnie zdefiniowany ciąg wielomianów ortogo- 
nalnych, bowiem w tym przypadku [a,b| = |-1,1] i waga jest jednost- 
kowa, p = 1. Mamy 


L „a 5 i 
Ly(z) = Ek dat © =1)5 
albo 
Lolz) = 1, L(x) = T, 
2k— 1 k—1 
L(£) = k cLy_1(1) — Lg-2(x), k >2. 


Wielomiany Czebyszewa Tę. 

Ten ciąg wielomianów poznaliśmy już w Rozdziale 7.2. Są one ortogo- 
nalne w przedziale [—1,1] z wagą p(z) = (1 — 2?) 12. Przypomnijmy, 
że Ty(x) = cos(k arccos(x)) albo 


T(x) = 1, BIE w 
T(x) 2xTp—ı(£) — Tk-2(£), k>2. 
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Wielomiany Hermite’a Hy. 

Jest to przykład wielomianów ortogonalnych na całej prostej rzeczywi- 
stej, a więc (a,b) = (—00, +00), przy czym waga p(x) = e77. Wielo- 
miany Hermite’a dane są wzorami 


albo 


H(z) = 1, Hı(x) = 2z, 
Hy(x) = 2xHy_,(x) — (2k — 2)Hk_o(z), kiz 2; 


Zauważmy, że we wszystkich tych przykładach kolejny wielomian 
ortogonalny zależał tylko od dwóch poprzednich. Nie jest to przypadek, 
ale ogólna reguła. Wielomiany ortogonalne spełniają bowiem formułę 
trójczłonową, zob. U. 10.1. 

Ciągi wielomianów ortogonalnych mają również szereg innych cie- 
kawych własności, co powoduje, że często używane są jako narzędzie 
do aproksymacji funkcji wielomianami, zob. U. 10.2. Teraz pokażemy 
tylko jedną bardzo ważną ich własność, a dokładniej, własność zer tych 
wielomianów, która będzie nam potrzebna do rozwiązania problemu 
kwadratur maksymalnego rzędu. 


Lemat 10.3 Wielomian ortogonalny pp ma dokładnie k pojedynczych 
i różnych żer w przedziale otwartym (a, b). 


Dowód Niech y,y2,...,ys będą wszystkimi punktami z przedziału 
(a,b), w których wielomian py zmienia znak. Gdyby lemat nie był praw- 
dziwy to mielibyśmy s < k. Wtedy, kładąc 


w(x) = (z — y)(t — ye) **-(0— ys), 
mamy (w, pk) = 0, bo w € Il, i s < k. Z drugiej strony, 
b 
(wp) = | w(a)pi(a)p(a)dz 7 0, 


bo funkcja podcałkowa jest albo stale dodatnia, albo stale ujemna (z 
wyjątkiem zbioru miary zero). Ta sprzeczność kończy dowód. 
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10.3 Kwadratury Gaussa 


Dla dowolnej n, niech zg, zj,... „z; będą (różnymi) zerami (n+ 1)-szego 
wielomianu p„+4 w ciągu wielomianów ortogonalnych w przestrzeni uni- 
tarnej L2 (a,b), tzn. 


Pne1(1) = (x — x0) (x — TÄ) (£ = £3). (10.1) 


Ponieważ x; leżą w przedziale (a,b), możemy mówić o kwadraturze 
interpolacyjnej opartej na tych węzłach. 


Definicja 10.3 Kwadraturę interpolacyjną QSS opartą na zerach wie- 
lomianu ortogonalnego pn+1 nazywamy kwadraturą Gaussa. 


Okazuje się, że właśnie kwadratury Gaussa mają najwyższy rząd. 
Dokładniej, mamy następujące twierdzenie. 


Twierdzenie 10.2 Kwadratura Gaussa QS? ma najwyższy rząd spo- 
śród wszystkich kwadratur opartych na węzłach o łącznej krotności n+1, 
oraz 

rz(Q$*) = 2n +2. 


Dowód Wobec Lematu 10.1(ii) wystarczy pokazać, że kwadratura QS5 
jest dokładna dla każdego wielomianu stopnia nie większego niż 2n + 1. 
Niech f € Ilanı. Niech w; € Il, będzie wielomianem interpolującym 
f w węzłach-zerach xġ,...,xž wielomianu p„+1. Jeśli deg f < n to 
oczywiście wę = f i QSS (f) = S,(wę) = S (f). Jeśli zaś 


n+1 < deg f < 2n+1, 


to f — wz jest wielomianem stopnia tego samego co f i zeruje się w zj, 
0<j <n. Stąd 


f(z) — wę(z) = (x — x)(x — xi) (x — zn)g(z), 
gdzie g jest wielomianem, 


deg g = deg f — (n +1) < (2n+1) — (n+1) =n. 
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Korzystając z (10.1) i faktu, że pn+ı jest prostopadły do IIn, ostatecznie 
otrzymujemy 


SA) - FO) = [ (A-vwOpOd 


II 
— 
a> 
R 
| 
R 
== 
w 
| 
8 
Z 
8 
3 
B 
m 
R 


co kończy dowód. 


Zajmiemy się teraz błędem kwadratur Gaussa. Pokażemy, że również 
ze względu na błąd mają one dobre własności. 


Twierdzenie 10.3 Jeśli f € C?"**([a,b|) to błąd kwadratury Gaussa 
QSS wyraża się wzorem 


AGIA n “E 
gdzie £ € [a,b]. Stąd, w szczególności, 
M |lPnll? 
S ZO = a 
jerżnt (ab) | (F) © (A) (2n + 2)! 


Dowód Niech wy € II, będzie wielomianem interpolującym f w zerach 
Tj wielomianu pn+1. Niech wy € Il,,, będzie z kolei wielomianem 
(Hermite'a) interpolującym f w dwukrotnych węzłach x, tzn. takim, 
że d;(z;) = f(x;) i w(x) = f'(x), 0 < j < n. Ponieważ rz(Q7*) = 
2n + 2 to QSF (Ŭŭ) = S (Ñr), a stąd i ze wzoru na błąd interpolacji 


Hermite'a mamy 


S„(f) — Gr (P) 


SAP) - Qu) = SF) — QE" (dy) 
SF) — Sly) = | (Fa) — By(0))pla) dz 


[e-r e- Paai opla) dz 


Il 
z 
I 
+ 
= 
B, 
D 
B 
=> 

8 
ox 
R 
3# 
B 

Q 

8 
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Ponieważ funkcja p, ,,(r)p(x) jest prawie wszędzie dodatnia, możemy 
teraz zastosować twierdzenie o wartości średniej, aby ostatecznie otrzy- 
mać 


| 72 


50) - QE) = N S Pa Opla)dz 


(2n+2) 
= E ol, 


co kończy dowód. 


Na końcu, zwrócimy jeszcze uwagę na inną, bardzo ważną własność 
kwadratur Gaussa, a mianowicie, że ich współczynniki są dodatnie. 
Rzeczywiście, zapisując 


OF) = Xafla) 


i podstawiając 


f;(a) = (z— zg)”: (2 — sj 4) (5— cj) ---(2— sy)” 


mamy, że fj € Ilp, i fj jest prawie wszędzie dodatnia. Stąd 
źSUSOF Ea) 


i az > 0, bo f,(x;) > 0. Przypomnijmy, że dodatniość współczynni- 
ków kwadratury ma duże znaczenie przy numerycznym ich obliczaniu, 
zwłaszcza gdy funkcja podcałkowa f ma stały znak, zob. Rozdział 2.5.2. 


Mimo niewątpliwych zalet kwadratur Gaussa, ich stosowalność ogra- 
niczają trudności w wyliczeniu pierwiastków wielomianów ortogonal- 
nych, gdy stopień wielomianu jest duży. Wyjątkiem są tutaj kwadra- 
tury interpolacyjne oparte na zerach wielomianów Czebyszewa, zob. U. 
10.5. 


Uwagi i uzupełnienia 


U. 10.1 Pokażemy teraz, że wielomiany ortogonalne £py |+>o w danej prze- 
strzeni £2„(a,b) spełniają następującą formulę trójczłonową. Załóżmy dla 
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uproszczenia, że współczynnik przy x* w wielomianie py jest dla każdego k 
jednością. Wtedy istnieją liczby 6% (dla k > 1) i y > 0 (dla k > 2) takie, że 


po(z) = 1, 
pi(z) = (z-A), (10.2) 
p(z) = (©— Bk)pPk-1(%) — Ypk-2(u), k22. 


Aby to pokazać zauważmy, że py, można dla k > 1 przedstawić w postaci 
rozwinięcia 


k—2 
Pr(2) = (© — ch1)pk-1(2) + 2 cjpj(z). 
j=0 
Mnożąc skalarnie obie strony tego równania przez p; dla 0 < s < k- 3, 
otrzymujemy 


0 = (pr(v),ps(t)) = ((£ = ck-1)Pk-1 (2), ps(£)) + cs(ps(t),ps(T)). 


Wobec tego, że (1—cy_1)ps(z) jest wielomianem stopnia mniejszego niż k—1, 
mamy 


((© — ck—1)pk—1(1),ps(1)) = (Pk-1(2), (£ — ck-1)ps(1)) = 0, 
a stąd c;(ps(1),ps(1)) = 0 i cs = 0. Możemy więc napisać, że dla k = 1 
mamy pı(x) = (z — (1), a dla k > 2, 
pk(z) = (x£ — Bk)pr-1(£) — YkPk-2(2), (10.3) 


gdzie Oy = Ck-1 Í Yk = Choo. Aby jawnie wyznaczyć 4 i yk, pomnożymy 
skalarnie obie strony równania (10.3) kolejno przez pķ—1 I pk-2. Otrzymujemy 


0 = (pk(t),Pk—1) 


((x — Bk)pk—1(7), pk-1(7)) — Yk(Pk-2(7),Pk-1(7)) 
(Tpk_1(©),Pk_1(1)) — Bk(pk-1(0),Pk-1(1)), 


8, = (zpk—1(©), pk-1(1)) 
(pk—1(7),pk-1(©)) 


0 = (pk(£), pk-2(x)) 


((x — Br)pk-1 (£), pk-2(7)) — Yk(Pk-2(7),Pk-2(7T)) 
(pk-1(7), CPR2(T)) — Yk(Pk-2(tT), pk-2(7)), 
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a stąd i z równości (pk1(1), xpk—-2(£)} = (pk-1(1),pk-1(0)), 


(pk—1(7), pk—1(0)) 
(pr-2(x), pk-2(£)} 


IR = 


Zauważmy, że z formuły trójczłonowej wynika w szczególności algorytm 
wyznaczenia ciągu wielomianów ortogonalnych. Wystarczy bowiem wyzna- 
czać kolejne współczynniki Gy i yk (obliczając odpowiednie iloczyny skalarne 
(cpę(c),pę(z)) i (pk(x),pę(zt))) i stosować wzór rekurencyjny (10.2). Do- 
dajmy jeszcze, że w obliczeniach numerycznych najlepiej jest przechowywać 
informację o ciągu (pz |k>o po prostu w postaci liczb Dk i Yk- 


U. 10.2 Wygodnie jest posłużyć się wielomianami ortogonalnymi w przy- 
padku, gdy chcemy znaleźć najlepszą aproksymację danej funkcji f wielomia- 
nem ustalonego stopnia n, i błąd mierzymy w normie przestrzeni £2 „(a,b). 
Rzeczywiście, jak wiadomo, najlepszą aproksymacją dla f w przestrzeni II, 
jest jej rzut prostopadły na II„. Ponieważ n + 1 początkowych wielomianów 
ortogonalnych pk, O < k < n, tworzy bazę w Iln, rzut ten wyraża się wzorem 


U. 10.3 Zachodzi następujące twierdzenie Łuzina. Niech przedział [a,b] bę- 
dzie skończony. Niech Qn(f) = Żj=0 a f (=f) będzie takim ciągiem kwa- 
dratur, że: 


(n) 


(i) wszystkie współczynniki a, sẹ dodatnie, 


(ii) rząd kwadratur Qn rośnie do nieskończoności gdy n + ». 


Wtedy ciąg Qn(f) zbiega do P f(x)dx dla każdej funkcji ciągłej f. 

Zauważmy, że twierdzenie to stosuje się do ciągu kwadratur Gaussa QSS, 
ale nie do ciągu kwadratur Newtona-Cotesa QT, ponieważ w tych ostatnich 
pojawiają się dla dużych n współczynniki ujemne. 

W praktyce, najczęściej stosuje się ciąg kwadratur interpolacyjnych opar- 
tych na zerach kolejnych wielomianów Czebyszewa, ponieważ zera te dane są 
jawnie i *zagęszczają się”, tzn. zera wielomianu Czebyszewa 74 są też zerami 
wielomianu 124. Powstające w ten sposób kwadratury noszą nazwę formuł 
Clanshow-Curtis'a. Są one w pewnym sensie uniwersalne, bowiem posiadają 
optymalną szybkość zbieżności n""+) w klasach Fxy((a,b|) dla dowolnych 
riM. 
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U. 10.4 Jeśli przedział całkowania jest skończony i waga jest jednostkowa 
to kwadratury Gaussa (a dokładniej kwadratury Legendre'a) można użyć 
do tworzenia kwadratur złożonych OR , gdzie k oznacza liczbę podprzedzia- 
łów. Łatwo widać, że dla f € PA, bl), błąd takiej kwadratury można 
oszacować przez 


> M(b — a)?” t3 ,142r+2 
| (F) Qrk = (2r + 2)! (z) , 
czyli jest on porównywalny do błędu *zwykłej” złożonej kwadratury interpo- 
lacyjnej. Jednak złożona kwadratura Legendre'a korzysta z dwa razy mniej- 
szej liczby węzłów. 


U. 10.5 Błąd złożonej kwadratury Legendre'a w klasie Fy;"'([a,b]) można 
podać dokładnie. Wystarczy wykorzystać wzory z Ćw. 10.4 i 10.5, aby otrzy- 


mac IS) QGS| c ( 1 ak 
max — Qrk| = — ; 
JEET (lat) £ n 


gdzie 
C = 


M2 (© + 2) --- (2r + 2) | 
(2r +3)! TO N 
i n = k(r + 1) jest ogólną liczbą węzłów na [a,b]. 


Cwiczenia 

Ćw. 10.1 Uzasadnić, że: 

(a) jeśli kwadratura jest dokładna dla dowolnych n+ 1 wielomianów tworzą- 
cych bazę w IIn, to jest ona rzędu co najmniej n + 1, oraz 


(b) jeśli kwadratura jest rzędu n + 1 to jest ona niedokładna dla każdego 
wielomianu stopnia dokładnie n + 1. 


Ćw. 10.2 Uzasadnić, że kwadratura prostokątów jest kwadraturą Legen- 
dre'a, natomiast żadna z kwadratur Newtona-Cotesa QT nie jest kwadra- 
turą Gaussa. 


Ćw. 10.3 Załóżmy, że dane są liczby Oy i yk definiujące ciąg wielomia- 
nów ortogonalnych przez formułę trójczłonową. Zaproponować ekonomiczny 
(tzn. o koszcie proporcjonalnym do n) algorytm obliczania wartości n-tego 
wielomianu ortogonalnego w danym punkcie z, wykorzystujący formułę trój- 
członową. 
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Ćw. 10.4 Niech zj, 0 < j < n, będą zerami (n + 1)-szego wielomianu 
Legendre’a. Pokazać, że 
1 
J œ- x0)}?( - 21)? (2 — 2n)? d2 
—1 


= zaj 1-2::-:n:(n+1) j 
2n +3 i(n+2)(n+3):::(2n + 2) 

Wskazówka. Wykorzystać fakt, że dla n-tego wielomianu Legendre'a mamy 

fi Zę(a)da = (2n +1). 


Ćw. 10.5 Niech 


u (f) = $ afe) 
j=0 


będzie kwadraturą interpolacyjną opartą na zerach (n+ 1)-ezego wielomianu 
Legendre'a. Niech —00 < a < b < wo. Pokazać, że wtedy kwadratura 
> b—a<- cj +1 
(M) = 7 2,/(a+ —6b—0)) 
j=0 


jest kwadraturą Gaussa opartą na n + 1 węzłach, dla całki na przedziale 
(a,b) z wagą jednostkową. Ponadto, jeśli f e C?"+2)([a,b]), to 


b z b—ay2n+3, 1-2:--(n+1) 42 fOrT2(£) 
_ AGSIN — | 
ILO (6) = ( 2 (eae aaa) (2n + 3)! 
Ćw. 10.6 Niech p=li-=wo<a<b< +w. Pokazać, że jeśli kwadratura 
Q oparta na n +1 węzłach jest rzędu r > n + 1, to odpowiadająca jej 
kwadratura złożona Qęą ma następującą własność. Jeśli f e C)([a,b]) to 
lyr-1(b— a)” 


IS) — GPL < (2) 


Ćw. 10.7 Niech cj, 0 < j < n, będą zerami (n + 1)-szego wielomianu 
ortogonalnego Legendre'a (tzn. na przedziale |-1,1] z wagą 1). Niech dla 
0<j<n, 
- | Eo E-z)E-2a) (0-2) 
Wj 
—1 (£j — z0) ++- (Lj — £j—1)(£j — Lj+1) +: (Ej — Zn) 

Pokazać, że jeśli f i g są wielomianami stopnia nie większego niż n, to ich 
iloczyn skalarny w £21(—1, 1), 


IF lea): 


r! 


0,9 = [| asd = Y wte): 
i > 
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Rozdział 11 


Iteracje dla równań 
liniowych 


Algorytmy rozwiązywania układów równań liniowych postaci 


> 


AŻ — b, 


gdzie A jest nieosobliwą macierzą rzeczywistą n x n, a b jest wektorem 
rzeczywistym w R”, które rozpatrywaliśmy w Rozdziałach 3, 4 i 5, na- 
leżą do grupy algorytmów dokładnych albo bezpośrednich. To znaczy, że 
po wykonaniu skończonej liczby dopuszczalnych operacji elementarnych 
dostajemy w arytmetyce idealnej dokładne rozwiązanie 


"<A 
W tym rozdziale zajmiemy się algorytmami iteracyjnymi rozwiązy- 
wania układów równań liniowych. Polegają one na tym, że, startując 
z pewnego przybliżenia początkowego o, konstruuje się ciąg kolejnych 
przybliżeń 
le AD oe k=L2ux 


które w granicy osiągają rozwiązanie dokładne, 


lim Zy = 7*. 


k—oo 
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11.1 Kiedy stosujemy iteracje? 


Jasne jest, że algorytmy iteracyjne stosujemy wtedy, gdy są one konku- 
rencyjne w stosunku do algorytmów bezpośrednich. Dlatego przekształ- 
cenia $y należy wybierać tak, aby kolejne przybliżenia można było ła- 
two obliczać i jednocześnie kolejne błędy ||ż, — x*'|| szybko zbiegały do 
zera. 

Zwykle zakłada się również, że dokładne rozwiązanie £” jest punk- 
tem stałym przekształcenia $©,(A, b; -). Wtedy kolejne błędy spełniają 
zależność 

Tp — T* = O(A, b; Zo) — D(A, b; 7”). 


Jeśli teraz ©,(A, b; -) są lipschitzowskie ze stałymi my < +oo, tzn. dla 
pewnej normy wektorowej ||- || mamy 


I*x(1) — PDI < mcy], YZ, y, 


to 
[Zr — *|| < mf — z|]. 


Warunek liMmk—>oœ my = O jest więc dostateczny na to, aby metoda była 
zbieżna dla dowolnego przybliżenia początkowego o, przy czym szyb- 
kość zbieżności zależy od tego, jak szybko mą maleją do zera. Dla więk- 
szości stosowanych metod $y jest funkcją liniową błędu początkowego, 
tzn. 

D,(A, b; Zo — 7*) = My(ży — £*), 


gdzie M; jest pewną macierzą. Wtedy jako mą można przyjąć normę 
tej macierzy, 
m, = |Mzl| = sup |Mżl| 


|-||=1 
Dla ilustracji, rozpatrzmy ogólną metodę iteracji prostej, w której 


Ze = Bi, 4 +€, (11.1) 


dla pewnej macierzy B wymiaru n x n i wektora € € R”. W tym 
przypadku 
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a stąd i z nierówności ||B*|| < ||B|[", mamy 
lz -= z < BIFEZ — z". 


Warunkiem dostatecznym zbieżności iteracji prostych jest więc ||B|| < 
1. Mówimy, że metoda jest zbieżna liniowo z ilorazem |B||. 


Przykład 11.1 Rozkładając macierz A = (a,,);,-, na sumę 
A=D+C, 


gdzie D jest macierzą diagonalną składającą się z wyrazów stojących 
na głównej przekątnej macierzy A, układ Ar = b jest równoważny 
układowi 

Di = -C + b, 
a stąd (o ile na przekątnej macierzy A nie mamy zera) otrzymujemy 
metodę iteracyjną 

ik = BaRa FG 
gdzie B =—D"C ic= D-1b, zwaną metodą Jacobiego. 

W metodzie Jacobiego warunek dostateczny zbieżności, ||B|| < 1, 
jest spełniony wtedy, gdy macierz A ma dominującą przekątną, tzn. 
gdy s 

2|a::| > X |a|, 1<i<n. (11.2) 
j=1 


Rzeczywiście, ponieważ wyraz (i, j) macierzy D"!C wynosi 0 dla i = j 
i Qij aii dla i F J, to 


n 
IDC] = max >, dazl/laz,] 
j=Lljzi 
n 


= max, „|assl/las;] KN 
J 


przy czym ostatnia nierówność wynika z (11.2). 


Inne przykłady iteracji prostych podane są w U. 11.3 i Ćw. 11.2. 


Zastanówmy się teraz nad złożonością metod iteracyjnych. Ponieważ 
możemy jedynie znaleźć pewne przybliżenie rozwiązania dokładnego £*, 
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przez złożoność metody będziemy rozumieli koszt kombinatoryczny ob- 
liczenia 1% z zadaną dokładnością e > 0. Dla uproszczenia założymy, że 
medoda jest zbieżna liniowo z ilorazem m. Zauważmy, że aby zreduko- 
wać błąd początkowy do s > 0, wystarczy wykonać k = k(e) iteracji, 
gdzie k spełnia 

m*||ży — £*|| < e, 


czyli 

p > los(1/E) — log(1/liżo — *"[|) 

- log(1/m) 
Liczba ta zależy więc w istotny sposób od błędu początkowego i (przede 
wszystkim) od stałej Lipschitza m, natomiast zależność od dokładności 
e i wymiaru n układu jest dużo mniej istotna. Zakładając, że koszt jed- 
nej iteracji wynosi c = c(n) (zwykle c(n) jest tym mniejszy, im mniejsza 
jest liczba niezerowych elementów macierzy A), złożoność metody jest 
proporcjonalna do 
log(1/e) 


0, log(1/m)” 


Stąd oczywisty wniosek, że metody iteracyjne warto stosować zamiast 
metod bezpośrednich w przypadku gdy 


e wymiar n układu Ar = b jest “duży”, oraz 


e macierz A układu jest *rozrzedzona”, tzn. ma stosunkowo nie- 
wielką liczbę elementów niezerowych, np. proporcjonalną do n. 


Układy o tych własnościach powstają często przy numerycznym roz- 
wiązywaniu równań różniczkowych cząstkowych. 

Zaletą metod iteracyjnych jest również ich prostota, przez co są one 
łatwe do zaprogramowania. 


11.2 Metoda Czebyszewa 


Zauważyliśmy, że ze względu na szybkość zbieżności metody iteracyjnej 
ważne jest, aby stałe lipschitzowskie mą odwzorowań $, malały jak 
najszybciej. Metoda Czebyszewa, którą przedstawimy w tym rozdziale, 
jest właśnie próbą minimalizacji tych stałych. 
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Przy opisie metody będziemy korzystać z następujących dwóch fak- 
tów, które sformułujemy jako lematy. 


Lemat 11.1 Jeśli macierz A jest symetryczna, to istnieje ortonor- 
malna w R” baza jej wektorów własnych Ej, 1 < j < n, ten. 


Ea 


a odpowiadające im wartości własne Aj, AG; = AEn sq rzeczywiste. 
Jeśli ponadto A jest dodatnio określona to A; są dodatnie. 


Lemat 11.2 Niech macierz A będzie symetryczna i niech 
Al > [A| > -e > JAn] 


będą jej wartościami własnymi. Wtedy 


> 
|Aa = sup FE 
Z0 zll2 


= fàl. 


Jeśli ponadto A jest nieosobliwa, tzn. |Aņn| > 0, to 


1 
A™|l = —. 


Lemat 11.1 jest znanym faktem z algebry liniowej, więc dowód po- 
miniemy. Dowód Lematu 11.2 podany jest w U. 11.1. 


Zakładamy, że macierz A wyjściowego układu równań jest syme- 
tryczna i dodatnio określona. 


A=A'>0Q, 
a jej wartości własne A, leżą w znanym przedziale [a,b], 


0)<AANL<''-<A=]|AJ>2 <b < oo. 


W metodzie Czebyszewa kolejne przybliżenia 1, konstruujemy tak, 
aby była spełniona równość 


(1A-1) = Wz(A)(fo—2"), (11.3) 
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gdzie W;y(4) jest macierzą, która wielomianowo zależy od A, 
k (k) 
j=0 
i ma możliwie małą normę. (Dalej będziemy również używać tego sa- 
mego symbolu Wy do oznaczenia wielomianu zmiennej rzeczywistej, 
tzn. dla t € R mamy W(t) = 2 5 0) Możemy więc napisać, że 
Tk = W/(A)% + (I — WZ(A))T*. 


Aby zapewnić konstruowalność przybliżenia %ęķ, musimy umieć obliczać 
(I —Wy(A))z*. Ponieważ 


k 
(I -WAE = -PaP AZ + (1-a$*)f* 
j=1 
Ka 4 3 
= -) aA + S 
j=0 
dla konstruowalności z, wystarczy założyć, że ag” = 1, co odpowiada 
warunkowi 
W,(0) = 1, Vk. (11.4) 


Z nierówności (11.3) mamy 


lE — ll < [Wz(A)|allfo — ©"||>. 


Idea metody polega teraz na wyborze wielomianów W/4 tak, aby miały 
one jak najmniejszą normę ||[W;4(4)||-. W tym celu zauważmy, że ma- 
cierz W;(4) jest symetryczna, a jej wartości własne wynoszą W;(A;), 
gdzie A; są wartościami własnymi macierzy 4, zob. Ćw. 11.1. Stąd 


[WA < max |Wz(t)| = |IW|lo(taż)- 


Rozwiązanie problemu minimalizacji ostatniego maksimum ze względu 
na wszystkie wielomiany W, € II, takie, że Wy(0) = 1, jest podane 
w Rozdziale 7 (zob. Ćw. 7.8). Przypomnijmy, że optymalny wielomian 
dany jest wzorem 


MOSO 


11.2. METODA CZEBYSZEWA 


gdzie Tę jest k-tym wielomianem Czebyszewa, a 


= 2t — (b + a) 


h(t) EE 


Przypomnijmy również, że wielomiany Czebyszewa spełniają formułę 
trójczłonową, Ty(t) = 1, Ty(t) = t, oraz Ty(t) = 2tTk-1(t)—Tk-2 dla k > 
2. Pozwala to na rekurencyjną konstrukcję kolejnych z4 w zależności od 


Tk=1 1 Tk-2. Dokładniej, dla pierwszego przybliżenia mamy 


2A—(b+a)I 
W:(A) — TKA) _ (=P NEI. 
i Tı (A(0)) (a) bpa 
a stąd 

2 z B DE wak- wyk 
11 = W, (A)zo + (I - W;(A))z = fo + ma To) 

= + 2 wd 

m war 


gdzie Tę = b — Az jest początkowym residuum. 


Dla k > 2 wykorzystamy formułę rekurencyjną dla wielomianów 


Czebyszewa. Oznaczając t; = T;(h(0)), j > 0, mamy 


T.(h(A)) _ 2h(A)TR_1(h(A)) — Tk-2(h(A)) 


w;(4) = = 


| 
N 
~œ 
w 
| L 
za~ 
— 
= 
— 


k E tk 
tk— ty 
= 257 h(dlózi=F) = E a 
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158 
tkı7 2 b+a tk-2 
tk C b—a ) E 2 tk (Zr-2 2 
Ee a e otk-10+4 z tk-2 z 
= DES E k-1 Ep sù k—1 ES k—2 
tkıb+a tk—2 
FAZ H T". 
( tk b—a tk )a 


Wykorzystując wynikającą z formuły trójczłonowej równość 


b—a 
b+a 


tk = —2 tk—-1 — tk-2 


otrzymujemy ostatecznie 
k = (1+ ak)fk-1 — Qklk-2 — OkTk-1, 


gdzie ap = tk—2/tk 1 Bk = (4tk-1)/(tz(b + a)). Algorytm wynikający z 
metody Czebyszewa można więc zapisać następująco. 


To := {dowolne}; fą := b — Afo; 
Tı := To + 2To/(b + a); M := b— Ad; 


p = —4/(b =, a); 
for J = 2,0, do 
begin 


tk = 2titk-1 — tk-2; 


ak := tk-2/tk; 
Tp := (1 + Gk)Tk=1 — Owik=2 + bltk-1/tr)fk-1; 


Te := b — Aż, 
end. 
Zauważmy, że opisana metoda Czebyszewa jest metodą dwupunktową, 
bowiem do konstrukcji kolejnego przybliżenia 1, wykorzystuje się dwa 
poprzednie przybliżenia Tę_1 i Tę». 
Zastanowimy się jeszcze nad szybkością zbieżności metody Czeby- 


szewa. Wobec 
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mamy 
[Zo — 7*|| 


nE 


b—a 


lZk — 7*ll2 < 


Aby oszacować T; '((b + a)/(b — a)), wykorzystamy jawną formułę na 
T(t) z Cw. 7.2. Oznaczając t = (b + a)/(b — a) mamy dla dużych k 


1 2 2 
TROA UFP- +0- (tve 
2 
= k 
(= | ei 
b+a | (ba)? 


2(b— a)* oz; 


a stąd w przybliżeniu 


k 
/b/a— 1 

lZk = zll < 2| = ] IZo z a 
/b/a+1 


Metoda Czebyszewa jest więc zbieżna co najmniej liniowo z ilorazem 
(4/B/a— 1)/(y/b/a + 1). Zauważmy jeszcze, że jeśli a i b są odpowiednio 
najmniejszą i największą wartością własną macierzy Á, to b = ||All2, 
1/a = || At], a stąd b/a = |A||>||A""||2 = cond(A) i metoda jest 
zbieżna z ilorazem 


cond(A) — 1 
cond(A) + 1 


Po raz kolejny widzimy tu ważną rolę uwarunkowania macierzy. Im 
uwarunkowanie jest gorsze tym gorsza zbieżność metody Czebyszewa. 
11.3 Metoda najszybszego spadku 


Tak jak poprzednio, będziemy zakładać, że macierz A układu jest sy- 
metryczna i dodatnio określona. Wtedy zachodzi następujący fakt. 
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Lemat 11.3 Rozwiązanie i* układu Ať = b jest jedynym wektorem 
minimalizującym formę kwadratową 


QZ) = SITAT g b, 
Dowód Zauważmy, że dla dowolnych r i ż mamy 
QE+2) = „(a | DUG EA 
= LAr — T6 + RGZU + IT AŻ) — Z b+ iaz 


= Q(Z) — z7”b— AT) + s7 AŻ 


przy czym w ostatniej równości wykorzystaliśmy symetrię macierzy A. 
Stąd dlar=f*i Z 0 dostajemy 


QE +2) = QUE") + ZAZ > QUE), 


co z kolei wynika z dodatniej określoności A. 


W metodzie najszybszego spadku konstruuje się kolejne przybliżenie 
Zy przesuwając się od 1y_, w kierunku przeciwnym do gradientu formy 
Q tak długo jak długo wartość tej formy maleje. (Przypomnijmy, że 
minus gradient wyznacza kierunek lokalnie najszybszego spadku.) Po- 
nieważ gradient Q w punkcie Zęy_1 wynosi 


ð á = 
( S0) adinei a 


OT; i=1 


metoda ma postać 
Tk = Tko1 + YkTk—1, (11.5) 


gdzie y, jest dobrane tak, aby zminimalizować Q(x,). Jak łatwo się 
przekonać, 


Q(T) = QlTk-1 + YTk-1) 
= 1 z x) > 
= Q(Zk-1) + 7 Tk-1ATk-1 + yri Ak = YTĘ_10 
JE. 5 T > 
= Ql£k1) + EVT 1 Ark-1 — Yrk=1Tk-1 


2 
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jest jako funkcja y parabolą, która przyjmuje minimum w 


T > 
Th—17k—1 


GER * 11.6 
Yk rI Afk ( ) 


Wzory (11.5) i (11.6) definiują już w pełni metodę najszybszego spadku. 
Pokażemy jeszcze zbieżność metody. W tym celu zauważmy, że 


T ERC 
Q(Tk) = Q(Tx_1) — UAE < Q(T-1). 


2 Tk-1ATk-1 


Ciąg Q(T) jest nierosnący i ograniczony z dołu przez Q(x"), a więc 
jest zbieżny. To wymusza 


T 2 y2 
s TĘ—17k—1 
lim GSUŻU = 0, 

k—co Tk_1ATk—1 
co z kolei oznacza, że ry = b — Ady zbiega do zera i w konsekwencji 


k—oo 


Uwagi i uzupełnienia 


U. 11.1 (Dowód Lematu 11.2) Jeśli (dh jest ortonormalną bazą wek- 


torów własnych macierzy A, to dla dowolnego wektora i = yi aë; mamy 


(AT, AT)2 = > o, A(ć;), 5 AG) 
2 


i=1 j=l 


EG 


- (DE aE) = Xabi Ea 
i=1 j=l 2 j=1 


Aii£]|2, 


a stąd ||Ax||e < |X,|||Z||2. Z drugiej strony, 
lAl = [éle = Allé = ll, 


a więc || Alļl2 = [A1]. 

Aby pokazać drugą część lematu zauważmy, że macierz A"" jest też sy- 
metryczna. Ma ona te same wektory własne co A, a wartości własne wynoszą 
1/A;, przy czym |1/An| > ++ 2 [1/24]. Stąd ||A7"||2 = [1/2n]. 
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U. 11.2 Pokazaliśmy, że warunkiem dostatecznym zbieżności metody ite- 
racji prostej tę = Biy_; + Z jest ||B|| < 1. Okazuje się, że można podać 
warunek konieczny i dostateczny. Metoda iteracji prostej jest zbieżna dla 
dowolnego przybliżenia początkowego 19 wtedy i tylko wtedy, gdy promień 
spektralny macierzy B jest nie większy od jedności, tzn. 


p(A) = max |A| < 1, 


gdzie maksimum wzięte jest po wszystkich wartościach własnych macierzy 
A. 


U. 11.3 Przedstawiając macierz A w postaci A = L + U, gdzie L jest ma- 
cierzą trójkątną dolną elementów występujących na i pod główną przekątną 
A, układ równań Ať = b jest równoważny układowi Li = -U + b. Stąd 
wynika metoda iteracji prostej 


Tk = -LİU + EB, 


zwana metodą Seidla. Można pokazać, że metoda ta jest zbieżna gdy macierz 
A jest symetryczna i dodatnio określona. 


U. 11.4 Pokazana metoda Czebyszewa dla macierzy A = AŤ > 0 o war- 
tościach własnych z przedziału [a,b], 0 < a < b < œ, jest metodą dwu- 
punktową. Rozpatruje się też jednopunktowe, cykliczne metody Czebyszewa 
postaci 
Tk = Tk-1 + QkTk-1, 

gdzie parametry a4 są wybierane cyklicznie, &j+m = a; dla wszystkich j i 
pewnej ustalonej m, oraz tak, aby zminimalizować normę macierzy przejścia 
w każdym cyklu. Dokładniej, mamy Zk — T* = (I — ay A)(Zk-1 — T“), a stąd 


Zm — T* = Wm(4A)(7o — £*), 


gdzie 

Wml A) = (T= amA)(I — am—14):::(I — 4). 
Wielomian W, (ż) jest więc stopnia m, a jego pierwiastki wynoszą t; = 1/a;, 
1 < j < m. Ponadto Wm(0) = 1. Jak już wiemy, norma macierzy |Wm(4)||2 
jest minimalna gdy Wm(t) = W(t) = Tm(h(t))/Tm(h(0)), gdzie h(t) = 
(2x — (a + b))/(b— a). Stąd parametry a;4m najlepiej jest wybrać tak, aby 
a; były odwrotnościami pierwiastków W, (t), czyli 


b+a b-a 2j— 1 
aj = 2% cos ( F 


=i 
»)) > 1<j<m. 
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U. 11.5 Metoda najszybszego spadku należy do rodziny metod gradiento- 
wych postaci 
Tk = Tk-1 + Qkľk—1, 


gdzie parametr ay, wybiera się tak, aby zminimalizować normę 


lZk- 2*|B = y (Er - ©)TB(G, — 7), 


a B jest pewną macierzą, B = BT > 0. Zwykle przyjmuje się B = 4?. Jak 
łatwo się przekonać, metodę najszybszego spadku otrzymujemy w przypadku 
p= 1, czyli B = A. 

Można pokazać, że dla metod gradientowych 


cond(4) — 1 


k 
cond( A) + z) [202a 


I-ze < ( 


Cwiczenia 
Ćw. 11.1 Niech A będzie macierzą symetryczną o wartościach własnych 


àj, 1 < j < n. Pokazać, że dla dowolnego wielomianu W(t) = Da ajti, 
macierz 


W(A) = X ajA 
j=l 


jest też symetryczna. Ponadto wektory własne A i W(A) są takie same, a 
wartości własne W (4) wynoszą W(A;) dla 1 < j <n. 


Ćw. 11.2 Niech A= AT > 0. Rozpatrzmy metodę Richardsona postaci 
ik = (I — aA)Tk—ı + ab. 


gdzie a dobieramy tak, aby zminimalizować normę macierzy ||I — aAl|o. 
Zauważyć, że jest to metoda iteracji prostej, ale równocześnie jednopunktowa 
metoda Czebyszewa z U. 11.4 przy m = 1. Pokazać, że optymalne a wyraża 


się wzorem 
2 


Qopt = EP 
gdzie Ay i X, są odpowiednio największą i najmniejszą wartością własną A. 
Ponadto 
Àl — An cond(A) — 1 


I — o A — = ; 
|Z = aorAl> N+An  cond(A)+1 


gdzie cond(A) = ||/A||2||A""||2. 
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Ćw. 11.3 Pokazać, że warunkiem koniecznym i dostatecznym zbieżności 
metody Richardsona z Cw. 11.2, dla dowolnego przybliżenia początkowego 
To, jest 


0<a< (11.7) 


2 
lAll2 
Ponadto dla a spełniającego (11.7) iloraz zbieżności wynosi 


> ad — 1 Qopt < 2/A1, 
[7 — aA|2 = 1—A0N, 0 < a < Qopt, 


gdzie &opt = 2/(M + An). 


Ćw. 11.4 Wykazać, że m-cykliczna jednopunktowa metoda Czebyszewa z 
U. 11.4 jest zbieżna liniowo z ilorazem równym (b — a)/(b + a) da m = 1 i 


dążącym do (vb — Va)/(vb + Va), gdy m — oo. 
Ćw. 11.5 Niech A= AT > 0 i niech 
Q(z) = z TAŻ + ab. 
Pokazać, że zbiór 
{ZE R”: Q(T) = const. | 


tworzy hiperelipsoidę w R” o środku A-!b i osiach w kierunkach wektorów 
własnych macierzy A o długościach proporcjonalnych do 1/,/X;, gdzie A; są 
wartościami własnymi A. 

Wskazówka. Wykorzystać Lemat 11.1. 


Rozdział 12 


Iteracje dla równań 
nieliniowych 


W Rozdziale 11 analizowaliśmy metody iteracyjne dla rozwiązywania 
równań liniowych Ax = b. Teraz zajmiemy się równaniami nielinio- 
wymi, w których zamiast operatora liniowego (macierzy) A mamy od- 
wzorowanie nieliniowe f : D — R”, gdzie D jest podzbiorem w R”. Dla 
uproszczenia będziemy zakładać, że b= Ù, gdyż w przeciwnym przy- 
padku możemy zastąpić f funkcją fı(Z) = f(z) — b. Bedziemy również 
zakładać, że f należy do pewnej znanej klasy F funkcji przekształcają- 
cych D w R” i mających co najmniej jedno zero w D. Zadanie polega 


więc na znalezieniu rozwiązania z* = £*(f) równania nieliniowego 
f(=0, dla  feF. (12.1) 


Jak zwykle w przypadku, gdy danymi są funkcje zakładamy, że jedyną 
informacją “a priori” o f jest, że f € F. Dodatkową informację mo- 
żemy zdobyć przez obliczanie wartości, a niekiedy i pochodnych f (o 
ile istnieją) w wybranych punktach. Zauważmy, że metody iteracyjne 
Czebyszewa i najszybszego spadku z Rozdziału 11 dla równań linio- 
wych korzystały właśnie z takiej (częściowej) informacji o macierzy A, 
chociaż formalnie zakładaliśmy, że dostępna jest pełna informacja o jej 
współczynnikach. Dla równań nieliniowych możliwość obliczenia jedy- 
nie informacji częściowej jest założeniem modelowym, a to wymusza 
już fakt, że każda metoda rozwiązywania równania (12.1) może dawać 
co najwyżej wynik przybliżony. 
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12.1  Bisekcja 


Metoda bisekcji, czyli połowienia, jest dość naturalną metodą obliczania 
zer skalarnych funkcji ciągłych określonych na danym przedziale [a,b] i 
zmieniających znak. Dokładniej, rozpatrzmy klasę funkcji 


F=flfecC(ab]) : f(a) <0<f6W)). (12.2) 


Oczywiście, każda funkcja f € F ma co najmniej jedno zero w [a,b]. 
Startując z przedziału [a, b|, w kolejnych krokach metody bisekcji obli- 
czamy informację o wartości f w środku przedziału, co pozwala nam, 
w zależności od znaku obliczonej wartości, zmniejszyć o połowę prze- 
dział, w którym na pewno znajduje się zero funkcji. Bisekcję realizuje 
następujący ciąg poleceń, po wykonaniu którego x jest przybliżeniem 
zera funkcji f z zadaną dokładnością e. 


BEZ 6 ZW: 

c := (a+b)/2; e:= (b—a)/2; 

while (e >e) do 

begin 
if (f(x) < 0) then zr := x else al := z; 
x = (zl + gr)/2; e:= e/2 

end. 


Z konstrukcji metody łatwo wynika, że po wykonaniu k iteracji (albo 
po obliczeniu k wartości funkcji) otrzymujemy z = 2%, które odległe jest 
od pewnego rozwiązania z* o co najwyżej 


|rk= z*| < (5) (>9). (12.3) 


Metoda bisekcji jest więc zbieżna liniowo z ilorazem 1/2. Choć ta zbież- 
ność nie jest imponująca, bisekcja ma kilka istotnych zalet. Oprócz jej 
prostoty, należy podkreślić fakt, że bisekcja jest w pewnym sensie uni- 
wersalna. Dla jej zbieżności wystarcza bowiem jedynie ciągłość funkcji. 
Poza tym możemy łatwo kontrolować jej błąd. Konsekwencją (12.3) jest 
bowiem następujący wniosek. 
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Wniosek 12.1 Dla znalezienia zera z dokładnością e > 0, wystarczy 
obliczyć w metodzie bisekcji 


wartości funkcji. 


Dodajmy jeszcze, że bisekcja minimalizuje błąd najgorszy w klasie 
F zdefiniowanej przez (12.2), wśród wszystkich algorytmów korzysta- 
jących z określonej liczby obliczeń wartości funkcji, zob U. 12.1. 


12.2 Iteracje proste 


Przedstawimy teraz metodę iteracji prostych dla rozwiązywania rów- 
nań nieliniowych, którą można traktować jako uogólnienie iteracji pro- 
stych dla równań liniowych. Najpierw równanie (12.1) przekształcamy 
do równania równoważnego (tzn. mającego te same rozwiązania) 


Z = o(ż). (12.4) 


Następnie, startując z pewnego przybliżenia początkowego xo, konstru- 
ujemy ciąg kolejnych przybliżeń 14, według wzoru 


Ty = (Tk), k > 1. 


Twierdzenie 12.1 Niech Do będzie domkniętym podzbiorem dziedziny 
D, 
Do = Do C D, 


w którym © jest odwzorowaniem zwężającym. To znaczy, o(Do) C Do, 
oraz istnieje stała 0 < L < 1 taka, że 


lo) — AG) < Lzy], Vz, y € Do. 
Wtedy równanie (12.4) ma dokładnie jedno rozwiązanie £*, oraz 


dla dowolnego przybliżenia początkowego čo € Do. 
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Dowód Wobec 


IF) = O(Fr2)|| < L ||Ek-1 = Zr-2ll 
-< Lt — zoll, 


[Zk — Zk- || 


IA 


dla k > s mamy 


k k 
lz- zl < 2, lg- ls 2 lz — ol 
j=s+1 j=s+1 


= Lo (Lg Lo + LEH- zoll < 


S 


1- L 


If — Toll. 


Ciąg fix |x jest więc ciągiem Couchy'ego. Stąd istnieje granica a = 
limk w Tk, która należy do Dg, wobec domkniętości tego zbioru. Po- 
nieważ lipschitzowskość 6 implikuje jej ciągłość, mamy też 


O eaa 


tzn. Œ jest punktem stałym odwzorowania 6. Dla jednoznaczności za- 
uważmy, że jeśliby istniał drugi, różny od æ, punkt stały 8 , to mieliby- 
śmy 5 5 p 

la Ell = ló(a) — ó(8)|I < LIl&- Al. 


Stąd 1 < L, co jest sprzeczne z założeniem, że $ jest zwężająca. © 


Z powyższych rozważań otrzymujemy natychmiastowy wniosek do- 
tyczący zbieżności iteracji prostych. 


Wniosek 12.2 Przy założeniach twierdzenia 12.1, metoda iteracji pro- 
stych jest zbieżna co najmniej liniowo z ilorażem L, tzn. 


| — £*|| < Z*|2,- "||. 
Dla ilustracji, rozpatrzmy natępujące proste równanie skalarne: 
z = cos(x), dla zLED=R. (12.5) 


W tym przypadku 6(x) = cos(x). Zauważamy, że w przedziale (0, 1] 
funkcja 9 jest zwężająca ze stałą 


R , gi 
L = max | cos (x)| = sin(1) < 1. 
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Stąd istnieje dokładnie jedno rozwiązanie równania (12.5) w przedziale 
[0, 1]. Rozwiązanie to może być aproksymowane z dowolnie małym błę- 
dem przy pomocy iteracji prostych, startując z dowolnego przybliżenia 
początkowego o € [0,1]. 

Zaletą iteracji prostych jest fakt, że zbieżność nie zależy od wymiaru 
n zadania, ale tylko od stałej Lipschitza L. Ma ona szczególne zastoso- 
wanie w przypadku, gdy funkcja 9 jest zwężająca na całym zbiorze D, 
tzn. Dy = D. Jeśli ponadto D ma skończoną średnicę diam(D), to dla 
osiągnięcia e-sproksymacji zera funkcji f wystarczy wykonać 


A _ Hog(lfo=£*[/e)1 _ flog(diam(D)/e) 
E OGAE JAA 127 4 


iteracji, niezależnie od zo. Metody zbieżne dla dowolnego przybliżenia 
początkowego, nazywamy zbieżnymi globalnie. 


12.3 Metody interpolacyjne 


Metody iterpolacyjne dla równań nieliniowych powstają podobnie jak 
kwadratury interpolacyjne dla całkowania. Najpierw daną funkcję in- 
terpolujemy wielomianem Lagrange'a albo Hermite'a ustalonego stop- 
nia, a następnie jako przybliżenie zera bierzemy zero tego wielomianu 
iterpolacyjnego. Proces taki można iterować (powtarzać) wymieniając 
za każdym razem jeden z punktów interpolacyjnych na punkt nowo 
obliczony. Zauważmy, że w ogólności takie metody muszą korzystać w 
kolejnych iteracjach nie tylko z ostatniego przybliżenia zę_4, ale również 
z kilku poprzednich, tzn. 


Tk — Ó(Tk-1, SE Ko EE 


Wymagają one również podania nie jednego, ale s przybliżeń począt- 
kowych £Tg,..., £s-1. 

W praktyce stosuje się jednak tylko interpolację wielomianami ni- 
skich stopni, np. wielomianami liniowymi lub najwyżej kwadratowymi, 
których zera można łatwo obliczyć korzystając z bezpośrednich formuł. 

Jak się przekonamy, metody interpolacyjne należą do grupy metod 
zbieżnych lokalnie. Znaczy to, że zbieżność ciągu (xy jy do zera danej 
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funkcji f jest zapewniona jedynie wtedy, gdy przybliżenia początkowe 
zostały wybrane dostatecznie blisko z*. 

W dalszych rozważaniach będziemy zakładać dla uproszczenia, że 
dziedzina D = R. 


W przypadku, gdy wielomianem interpolacyjnym jest wielomian li- 
niowy Hermite a, mówimy o metodzie Newtona, albo stycznych. Star- 
tując z pewnego przybliżenia początkowego £o, w kolejnych krokach 
metody, k-te przybliżenie z, jest punktem przecięcia stycznej do wy- 
kresu f w punkcie zy_1. Ponieważ równanie stycznej wynosi y(x) = 
fck) + F'(ck-1)(7 — £k-1), otrzymujemy wzór 


AE 20 
f' (Ga) 


Oczywiście, aby metoda Newtona była dobrze zdefiniowana, musimy 
założyć, że f'(cę_1) istnieje i nie jest zerem. 


Tk = Tk-1 — 


Zauważmy, że metodę Newtona można traktować jako szczególny 
przypadek iteracji prostych, gdzie 


$(z) = z — 


Widać też, że nie jest ona zbieżna globalnie. Nawet jeśli pochodna w 
zpk- Się nie zeruje, ciąg {£p} może nie zbiegać do zera funkcji f, zob. 
Ćw. 12.2. Okazuje się jednak, że jeśli wystartujemy dostatecznie blisko 
rozwiązania z*, to metoda Newtona jest zbieżna. Dokładniej, załóżmy 
najpierw, że f(x") = 0 oraz 


f(a*) # 0. 
Ponadto załóżmy, że f jest dwukrotnie różniczkowalna w sposób ciągły, 
f € C?(D). Rozwijając $ w szereg Taylora w punkcie z* otrzymujemy 
gx—2* = $(ai1) — ola”) = (zk1—«')$'(1*) + (zk1—2*)*$'(6k)/2, 


gdzie min(r*, £ķ-1) < k < max(r*,zy_1). Wobec tego, że $'(x*) = 


FEF" E/E)? =016'(4.) = F"E) F (Ek), mamy 


PME) 
2f'(6k) 


Tk — T* = (£k-1 — T“) (12.6) 
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Zdefiniujmy liczbę 


eO 
F'(2) 

Oczywiście R; > 0. Dla zy_, spełniającego |1y1—1x'| < R < Ry, 
mamy z równości (12.6) 


R; = sup sup 


r>0 fx:|r—z*|<r) 


KKJESICZEZAE 


gdzie q < 1 i q zależy tylko od R. 
Niech teraz x* będzie zerem m-krotnym, 


F) = fla) = = fM(')=07 Ma"), 


gdzie m > 2, oraz niech f będzie m-krotnie różniczkowalna w sposób 
ciągły. Wtedy 


Tk —T* = (£k-1— T*) — 


= (m1 -2)(1- >), (12.7) 


o ile zı jest “blisko” x*. 

Metoda Newtona jest więc zbieżna lokalnie. Z (12.6) i (12.7) można 
też wywnioskować, jaki jest charakter zbieżności metody Newtona. Dla 
zera jednokrotnego z* oraz f"(x*) #0 mamy bowiem 


3 ada 

2f'(z*) 
Mówimy, że zbieżność jest kwadratowa. Jeśli zaś f”(x*) = 0 to zbieżnośc 
jest nawet szybsza. Z kolei dla zera m-krotnego zbieżność jest liniowa 
z ilorazem (1 — 1/m). 


(cz — z”) m (c — Ep1) 


Metoda Newtona jest pierwszą poznaną tutaj metodą iteracyjną, 
która jest (dla zer jednokrotnych) zbieżna szybciej niż liniowo. Dla ta- 
kich metod wprowadza się pojęcie wykładnika zbieżności, który jest 
zdefiniowany następująco. 
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Powiemy, że metoda iteracyjna 6 jest w klasie funkcji F rzędu co 
najmniej p > 1, gdy spełniony jest następujący warunek. Niech f € 
F i f(x*) = 0. Wtedy istnieje stała C < o taka, że dla dowolnych 
przybliżeń początkowych £tg,...,%s_1 dostatecznie bliskich z*, kolejne 
przybliżenia zy = $(Tk_1,...,ty_s) generowane tą metodą spełniają 


|z — z*| < C|zpq—vx*|?. 
Ponadto, jeśli p = 1 to dodatkowo żąda się, aby Č < 1. 


Definicja 12.1 Wykładnikiem zbieżności metody iteracyjnej ġ w klasie 
F nazywamy liczbę p* zdefiniowaną równością 


p* = supfp> 1: 6 jest rzędu co najmniej p }. 


Możemy teraz sformułować następujące twierdzenie, które natych- 
miast wynika z poprzednich rozważań. 


Twierdzenie 12.2 Wykładnik zbieżności metody Newtona (stycznych) 
wynosi p* = 2 w klasie funkcji o żerach jednokrotnych, oraz p* = 1 w 
klasie funkcji o zerach wielokrotnych. 


Inną znaną i często używaną metodą interpolacyjną jest metoda 
siecznych, w której stosuje się liniowy wielomian interpolacyjny La- 
grange'a. Metoda ta wykorzystuje więc do konstrukcji x, przybliżenia 
£k—-1 1 Tk-2. Musimy również wybrać dwa różne punkty startowe £o i 
11. Ponieważ prosta interpolująca f w £k—1 1 £k-2 Ma Wzór 


otrzymujemy 


Wk-1 7 Uk-2 


fo) — f (£k-2) 


Zauważmy, że jeśli £ķ—1 I £ķ—2 są blisko siebie, to zy jest podobny do 
tego z metody Newtona, bowiem wtedy iloraz różnicowy (f(Tk-1) — 


Tk = Tk-1 7 


f(Gk-1). 
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f(£k-2))/(£k-1 — £k-2) dobrze przybliża pochodną f'(zy-1). Nie wy- 
starcza to jednak, aby osiągnąć zbieżność z wykładnikiem 2. Dokład- 
niej, można pokazać, że wykładnik zbieżności metody siecznych dla zer 
jednokrotnych wynosi p* = (1 + v5)/2 = 1.618..., zob. U. 12.3. 

Niewątpliwą zaletą metody siecznych jest jednak to, że nie wymaga 
ona obliczania pochodnej funkcji (co w praktyce jest często bardzo 
trudne, a niekiedy nawet niemożliwe), a tylko jej wartości. 


12.4 Obliczanie zer wielomianów 


Obliczanie zer w przypadku, gdy funkcja f jest wielomianem jest wła- 
ściwie osobnym zadaniem-my tylko dotkniemy tego tematu. 

Załóżmy, że wielomian w : R > R jest dany poprzez jego współ- 
czynniki a; w bazie potęgowej, 


n 
w(x) = X ajz’ 
j=0 


(an 4 0), i zastosujmy metodę Newtona dla znalezienia największego 
zera tego wielomianu. Jedna iteracja metody Newtona jest wtedy sto- 
sunkowo prosta, gdyż obliczenie wartości i pochodnej wielomianu w ko- 
lejnym pukcie € = zy_+ polega po prostu na dwukrotnym zastosowaniu 
schematu Hornera. Rzeczywiście, przypomnijmy, że schemat Hornera 
oblicza nie tylko wartość w(£), ale również współczynniki w bazie potę- 
gowej ilorazu wy(z) przy dzieleniu w(x) przez jednomian (x — 6), zob. 
Ćw. 6.2. Zauważmy dalej, że równość 


w(x) = (z — 6)w,(x) + w(ś) 


implikuje w'($) = w (£). Stąd schemat obliczania jednocześnie w(ć) i 
w'($) wygląda następująco. 


w0 := wl := a; 
for j := n—1 downto 1 do 
begin 
w0 := w0 * Ẹ + aj; 
wl := wl * E + w0 
end; 


w0 := w0 + av. 
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Po wykonaniu tych poleceń mamy w0 = w(ć) i w1 = w'(6). 


Twierdzenie 12.3 Załóżmy, że wielomian w(x) = X;-, aja) ma zera 
rzeczywiste 
E d 22 d én. 


Jeśli tylko przybliżenie początkowe £o > 41, to ciąg zy generowany me- 
todą Newtona jest ściśle malejący i zbieżny do €,. 


Dowód  Wielomian w można przedstawić w postaci 
w(x) = (x — ś1)w: (u), 
gdzie w;(1) = an(x — 62) :*: (1 — En). Wtedy 
w (z) = w(x) + (r — &1)wi (z), 
a stąd 
(12.8) 
Wi (£p) | , 


Wi(Tk1) + (Tk—1 — E1 )Ww1 (£k-1) 


Żk=Gr = (t17) Wara) 


(£k-1 = 4) ( z 


Zauważmy teraz, że wielomiany w(x) i w{(x) mają stały znak dla 
c > 4, który równy jest znakowi a,. Rzeczywiście, ponieważ w;(x) 
ma wszystkie n — 1 pierwiastków w przedziale (—00, £], to w4(xz) ma 
również wszystkie n — 2 pierwiastki w tym samym przedziale. Stąd oba 
wielomiany mają ten sam znak na (4,, +oo) równy znakowi współczyn- 
nika przy najwyższych potęgach tych wielomianów. Wspłczynnik ten 
jest w obu przypadkach równy an. 

Z rozważań tych wynika, że mnożnik przy (zy-1 — 4,) w prawej 
stronie równania (12.8) jest dla zķ-ı > 4, dodatni i mniejszy od 1, 
a stąd, że ciąg zy jest malejący i ograniczony z dołu przez é. Niech 
g = lim, wozy > 4. Biorąc granicę obu stron (12.8) przy k * œœ 
otrzymujemy 


wi(g) | | 


g—& = (GG) h — wi(g) + (g — G1)wi(g) 
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co wymusza g = é, kończąc dowód. © 


Dodajmy, że jeśli zero 6, jest jednokrotne, é > £2, to mamy zbież- 
ność kwadratową. Należy jednak podkreślić, że zbieżność kwadratową 
uzyskujemy dopiero pod koniec procesu iteracyjnego. Mamy bowiem 


! w(x) SE: 
-n "MERC Ib 


Jeśli więc wystartujemy z £o >> 4q to początkowo zaobserwujemy zbież- 
ność tylko liniową z ilorazem 1— 1/n. W miarę zbliżania się do 4, zbież- 
ność będzie rosła, aż osiągnie w końcu poziom zbieżności kwadratowej, 


wi (51) 
Tp — Gy FE (Tk-1 — 41)” , 
D wi (61) 
Dla zera wielokrotnego, 6, = *:: = Êm > śŚmi11, m < n, zbieżność pozo- 


stanie jednak cały czas na poziomie liniowym, przy czym iloraz zbież- 
ności wzrośnie do I — 1/m. W końcu, jeśli é jest zerem o maksymalnej 
krotności n, to łatwo zauważyć, że 


Tk—& = (1 = HN (era — &1). 


Pozostaje jeszcze problem znalezienia przybliżenia 19 > 44. Można 
to zrobić na kilka sposobów. Na przykład, zauważmy, że dla 


r> M= max (1, lasl) 
j=0 lanl 
mamy 
w(x) a S layl 7 1 © [ajl 
> — „SS .. ZEE LI 
eR (1 > e > Mż; mb SÜ: 


Stąd w(x) ma stały znak na (M, +00), a to oznacza, że 4, < M. Jako 
zo można więc wziąć dowolną liczbę większą od M. 
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Uwagi i uzupełnienia 


U. 12.1 Metoda bisekcji jest optymalna w następującym sensie. Niech A: 
F — R będzie dowolną metodą (algorytmem) aproksymującą zero zx*(f) 
funkcji f z klasy F zdefiniowanej w (12.2), korzystającą jedynie z obliczeń 
(informacji o) f w k punktach. Wtedy dla dowolnego y > 0 istnieje funkcja 
fy € F mająca tylko jedno zero x*(f,) w [a,b] i taka, że 


A). > (5) (57) +7. 


Co więcej, można pokazać, że fakt ten zachodzi też w węższej klasie F3 
funkcji f € F, które są dowolnie wiele razy różniczkowalne. Oczywiście, 
nie wyklucza to istnienia metod iteracyjnych (takich jak metoda Newtona), 
które dla f € Fı są zbieżne szybciej niż liniowo. 


U. 12.2 Wszystkie metody iteracyjne opisane w tym rozdziale należą do 
grupy metod adaptacyjnych. Są to metody, które obliczają kolejną wartość 
(ew. pochodną) funkcji f w punkcie zy, którego wybór istotnie zależy od po- 
przednio obliczonych wartości f, tzn. zy = ay(f(%0),..., f(ck-1)). Okazuje 
się, że metody nieadaptacyjne są dużo gorsze przy rozwiązywaniu równań 
nieliniowych. Na przykład, najlepsza metoda korzystająca z k wartości funk- 
cji wybranych nieadaptacyjnie ma w klasie (12.2) błąd najgorszy równy tylko 
(b — a)/(2k + 2). 

Zauważmy jeszcze, że dla zadania całkowania z Rozdziałów 9 i 10, opty- 
malne metody okazały się nieadaptacyjne. 


U. 12.3 Uzasadnimy teraz, że dla zer jednokrotnych wykładnik zbieżności 
metody siecznych wynosi p* = (1 + v5)/2. 
Niech f € C?(D). Wykorzystując ilorazy różnicowe otrzymujemy 


(Ga) 
f(2k-2, £k-1) 
f(£k-1,2*) ) 
IC Z 0) 
f(Ek-2, £k-1, 2“) 
f(2k-2, £k—1) 


zk" = (£k-1— T*) 


= (zk-1— T“) (1 


= (gp 56 (ij — T*) 


(gk — ©')(wk-2— W”) 
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Stąd widać, że dla £ķ—2 i 1ę_1 dostatecznie bliskich c”, ciąg zę zbiega do 1”. 
Aby znaleźć wykładnik zbieżności, załóżmy bez zmniejszenia ogólności, że 
f"'(x*) 4 0 (w przeciwnym przypadku zbieżność jest jeszcze szybsza). Wtedy 


Tk — L“ m C(ay-1 — z*)(zky-2 — x*), (12.9) 


gdzie C = f"(x*)/(2f'(x*)). Stąd, po podstawieniu z; = ln(|C(z; — x*)|), 
dostajemy 

Zk R Zk—-1 +F Żk—2. 
Ostatnie równanie jest równaniem różnicowym, którego rozwiązaniem jest 
ży m q, gdzie q = (1 + V5)/2 jest dominującym pierwiastkiem równania 
2? =z +1. Ostatecznie więc |xy — z*| £ |CTtet'|, czyli 


ez = a*l = OTE (10e) = [cepa = at. 


U. 12.4 Dla znalezienia zera funkcji f : R” — R” można zastosować wielo- 
wymiarową metodę Newtona. Jest ona zdefiniowana następującym wzorem 
rekurencyjnym: 

Zr = dwa — (F(Żxa)) "/(ik-1), 


gdzie f'(6) jest macierzą” Jacobianem funkcji f wziętym w punkcie £, tzn. 


AT) 5 (6) 
f(Ó = : i 
G(E) |. aë 
f1(6) Ś1 
JOE ls s 
PO En 


W praktyce metodę tą realizuje się rozwiązując w każdym kroku iteracyjnym 
układ równań liniowych n x n. Dokładniej, 


k = dk-1 — Uk—1, 
gdzie yk_1 jest rozwiązaniem układu 
f Gro)dk1 = F(Zk-1). 


Tak jak w przypadku jednowymiarowym, wielowymiarowa metoda Newtona 
jest zbieżna lokalnie do zera £* z wykładnikiem 2, o ile Jacobian w punkcie 
z* jest nieosobliwy, tzn. det(f'(£*)) £ 0. 
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Cwiczenia 


Ćw. 12.1 Uzasadnić, że metoda iteracji prostych zy = cos(cę_1) jest zbie- 
żna do jedynego rozwiązania równnania z = cos(x) dla dowolnego przybli- 
żenia początkowego xo € R. 


Ćw. 12.2 Podać przykład funkcji f i przybliżenia początkowego xo takich, 
że ciąg generowany metodą Newtona oscyluje, tzn. 125 = £o 1 Tos+1 = £1, 
dla wszystkich s. 


Ćw. 12.3 Jeśli w modelu obliczeniowym liczenie pierwiastka kwadratowego 
va (a > 0) nie jest dopuszczalne, to możemy go przybliżyć stosując wzór 
rekurencyjny 
(214 z) 
Tk = = (Tki F , 
2 Tk—1 


który powstaje przez zastosowanie metody Newtona do równania f(x) = 


x?—a = 0. Pokazać, że ciąg zy jest zbieżny do v/a dla dowolnego przybliżenia 


początkowego xo > 0. Scharakteryzować szybkość zbieżności. 


Ćw. 12.4 Zaproponować metodę iteracyjną obliczania 1 /a dla a > 0 nie 
używającą dzielenia. Jak wybrać przybliżenie początkowe aby metoda była 
zbieżna? Jaki jest wykładnik zbieżności? 

Wskazówka. Zastosować metodę Newtona do funkcji f(x) = 1/z — a. 


Ćw. 12.5 Rozpatrzmy metodę iteracyjną 
z k-17 4 
e S (aka) — Fla) 
gdzie a jest pewnym parametrem. Pokazać, że metoda ta jest zbieżna liniowo, 
o ile a jest dostatecznie blisko zera z*. 


f(£k-1), 


Ćw. 12.6 Pokazać, że metoda Steffensena 


fla) 
JOKE f(£k-1)) — f (£k-1) 


jest zbieżna lokalnie z wykładnikiem 2. 


Tk = Tk-1 


Ćw. 12.7 Rozpatrzmy metodę iteracyjną daną wzorem 


f(ck1) 

Pi: 

Pokazać, że jeśli f ma m-krotne zero z* to metoda ta jest lokalnie zbieżna 
do z* z wykładnikiem 2. 


Tk = Tk-1 M 


